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1  IDENTIFICATION  AND  SIGNIFICANCE  OF 
THE  PROBLEM 

Robust  control  is  a  method  of  maximally  enhancing  stability  and  performance  of  a 
dynamic  system  in  the  presence  of  modeling  and  environmental  uncertainties.  It  is 
natural  to  consider  using  robust  control  methods  to  improve  the  dynamic  behavior  of 
Large  Space  Structures  (LSS). 

A  key  aspect  of  robust  control  is  adaptability.  Adaptability  is  enhanced  when  there 
exists  a  mechanism  for  updating  the  control  dynamic  model  either  on  a  periodic  and/or 
an  as-needed  basis.  The  need  for  internal,  automated  model  updates  becomes  more 
acute  when  either  the  dynamic  plant  is  difficult  to  model  or  when  its  operating  environ¬ 
ment  is  either  poorly  defined  or  conducive  to  causing  changes  in  the  dynamics  of  the 
physical  plant. 

These  conditions  apply  to  the  LSS  in  orbit.  In  the  first  place,  models  of  space 
structures  which  sufficiently  describe  the  dynamics  are  quite  often  of  unrealistically 
high  system  order.  Researchers  are  now  analyzing  ways  to  meaningfully  utilize  models 
of  order  10,000.  A  great  deal  of  research  has  been  devoted  to  model  approximation 
schemes,  and  it  will  likely  remain  a  major  activity  well  into  the  future,  as  many  schemes 
are  application  dependent.  A  related  problem  is  that  analytic  models  are  hard  to 
verify  on  hardware  on  the  1-g  Earth.  It  appears  that  a  scheme  consisting  of  online 
identification  plus  a  controller  which  adapts  optimally  to  updated  dynamic  information 
is  a  meaningful  basis  for  highly  robust,  active  LSS  control. 

The  crux  of  the  problem  is  that  LSS  dynamics  are  essentially  infinite-dimensional, 
best  described  by  the  mathematics  of  Hilbert  spaces  and  partial  differential  equations. 
In  addition,  classic  and/or  reliable  control  synthesis  methods  generally  require  finite¬ 
dimensional,  linear  models  (Balas,  1982).  Further,  the  realities  of  processing  throughput 
requirements,  not  to  mention  software  and  hardware  issues  of  reliability,  redundancy 
management,  etc.,  create  a  strong  need  for  using  control  design  models1  whose  system 
order  is  very  judiciously  selected. 

Finding  an  appropriate  scheme  (or  family  of  schemes)  for  LSS  control  is  becom¬ 
ing  an  increasingly  critical  issue  as  national  policy  moves  toward  deployment  of  these 
structures  in  space  for  many  purposes,  including  scientific  research,  communication, 
manufacturing,  manned  life  support,  and  defense.  Defense  and  scientific  requirements 
in  particular  place  very  tight  performance  requirements  on  such  metrics  as  pointing 
accuracy.  Pointing  accuracy  is  even  more  difficult  to  achieve  due  to  one  effect  of  the 
economics  of  space  operations:  lighter,  more  flexible  and  larger  structures.  Effective 
control  of  slewing,  pointing  and  other  LSS  attitude  maneuvers  clearly  implies,  then, 
control  of  structural  flexibility. 

Typical  LSS  performance  metrics  -  for  example,  beam  or  telescope  line-of-sight 


1For  the  purposes  and  scope  of  this  report,  three  classes  of  LSS  models  are  referred  to  in  the  text:  (l) 
“truth”  models  -  reasonably  accurate  simulation  models  of  an  actual  LSS  structure;  (2)  control  design 
models  -  of  lesser  order,  used  in  control  synthesis;  (3)  on-board  models  -  part  of  the  flight  software, 
typically  in  the  controller,  estimator  and/or  identifier  subsystems.  Often,  but  not  necessarily,  on-board 
models  are  equivalent  to  control  design  models.  Model  types  (2)  and  (3)  are  usually  linear,  and  (1)  may 
be  nonlinear. 


(LOS)  accuracy  (ranging  near  1  /zrad),  modal  damping,  surface  shape  errors  (antenna 
applications),  etc.  -  can  be  highly  sensitive  to  changes  in  the  LSS  model,  as  well  as  to 
sensor  measurement  and  to  sensor  and  actuator  location.  Another  performance  issue 
not  usually  very  critical  in  the  case  of  design  exercises  operating  on  simpler  dynamic 
models  (e.g.,  rigid  body)  concerns  the  dynamic  interactions  of  LSS  sensors  and  actua¬ 
tors  with  the  main  structure.  These  interactions  arise  not  only  from  the  nature  of  these 
devices,  but  also  from  their  location  on  the  LSS. 

The  point  of  this  discussion  is  that  the  LSS  control  design  problem,  of  its  nature, 
has  extra  dimensions  of  interrelationship  among  the  modeling,  dynamics  and  control 
synthesis  disciplines,  than  are  found  in  most  other  control  design  problems. 

Another  implication  of  the  discussions  above  is  that  “simple”  (ie,  low  order)  dy¬ 
namic  models  are  traditionally  nearly  always  appropriate  for  initiating  a  design  effort 
or  establishing  proof  of  concept  to  a  control  design  idea.  This  approach  is  undertaken 
at  much  greater  risk  in  the  case  of  LSS  control  design.  The  so-called  spill-over  effects 
which  arise  from  the  dynamic  interaction  of  unmodeled  with  modeled  modes  pressure 
the  designer  to  use  a  higher  order  design  model,  while  design  algorithmic,  software  de¬ 
velopment,  processing  hardware,  system  throughput  and  reliability  issues,  etc.,  exert 
strong,  opposite  pressures.  This  is  a  classic  tradeoff,  but  it  is  very  amplified  in  the  LSS 
problem.  Reasonable  LSS  dynamic  models  usually  have  far  too  many  modes  to  admit 
effective  control  design  and  implementation. 

Briefly,  modeling,  dynamics  and  control  are  highly  interrelated  design  disciplines. 
Control  design  is  driven  primarily  by  spacecraft  dynamic  characteristics,  performance 
requirements,  type  and  degree  of  external  and  internal  disturbances,  and  type,  number 
and  availability  of  sensors  and  actuators.  The  elastic  or  bending  modes  are  theoretically 
infinite  in  number,  and  typically  have  low  natural  damping.  Performance  criteria  must 
be  well  defined,  and  all  disturbances  must  also  be  well  defined  and  modeled,  in  order  to 
have  a  well-posed  control  problem.  Limitations  on  control  bandwidth  dictate  a  control 
design  based  on  a  reduced  order  model  (ROM:  the  so-called  “control  design”  model), 
which  contains  only  selected  low  frequency  modes.  A  critical  design  issue  to  be  addressed 
is  efficient  identification  and  selection  of  these  ROM  modes. 

The  key  low  frequency  modes  usually  require  the  most  control  effort  due  to  their 
larger  settling  time,  and  to  their  larger  effect  on  performance  degradation.  However, 
the  design  method  ought  not  allow  the  more  poorly  known  high  frequency  modes  to  go 
unstable. 

The  Phase  I  effort  of  this  work  has  explored  the  feasibility  of  applying  robust  identi¬ 
fication  and  control  techniques  to  the  LSS  control  problem.  A  more  detailed  description 
of  these  techniques,  and  also  of  their  application  to  some  relatively  small  albeit  repre¬ 
sentative  LSS  models,  is  provided  in  later  sections.  In  the  following,  refer  to  Figure  1. 

The  identification  method  applies  a  canonical  variate  analysis  (CVA)  technique  to 
selected  input  and  measurement  signals  of  the  process  to  be  identified,  and  thereby  ex¬ 
tracts  model  parameters.  Model  order  is  statistically  determined  using  built-in  criteria. 
By  using  such  computationally  stable  algorithms  as  singular  value  decomposition,  CVA 
identification  is  well-suited  to  LSS  applications.  In  fact,  its  major  applications  success 
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Figure  1:  LSSICS  Overview,  Showing  Functional  Relationship  of  Modeling,  Dynamics 
and  Control 


to  date  is  in  the  very  related  area  of  flutter  suppression  (Larimore  and  Mehra,  1985), 
in  which  real  time  identifications  were  performed  successfully  on  a  full  model,  1/4  scale 
F-16  with  wing  stores,  in  a  wind  tunnel  at  NASA/Langley.  This  was  accomplished  in 
a  complete  “hands  off"  mode.  CVA  identification  is  discussed  in  more  detail  in  Section 
2.2. 

The  control  approach,  called  Model  Predictive  Control  (MPC),  is  a  multi- input/multi¬ 
output  (MIMO)  output  controller,  capable  of  driving  several  outputs  simultaneously  to 
their  individual  (input  commanded)  setpoints,  or  of  tracking  an  input  reference  trajec¬ 
tory.  The  control  needed  to  achieve  the  desired  output,  or  response,  values  is  obtained 
by  comparing  the  difference  between  a  “desired”  output  reference  trajectory  and  the 
predicted  no-input  trajectory,  against  the  predicted  control  response  trajectory.  The 
control  equation  can  be  converted  into  an  overdet ermine d  system,  and  solved  for  the 
control  inputs  using  a  weighted  minimization  criterion.  These  weights  become  con¬ 
trol  design  parameters,  in  a  manner  similar  to  weighted  linear  quadratic  (LQ)  design. 
Other  design  parameters  relate  to  shaping  the  output  trajectory  and  to  the  length  of 
the  control  prediction  window.  See  Section  2.3  for  more  detail. 

Most  standard  control  design  methods  have  difficulty  in  the  LSS  arena  because  they 
work  best  on  single-input-single-output  (SISO)  systems,  or  if  MIMO,  have  difficulty 
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in  dealing  with  spill-over  dynamics,  failed  actuators,  changing  dynamics,  and  control 
energy  limits.  The  model  predictive  control  scheme  proposed  here  as  a  viable,  (logically) 
flexible,  robust  MIMO  method  is  able  to  address  all  of  these  concerns  as  an  implicit 
feature  of  its  design  process.  MPC  by  itself  is  robust;  it  is  therefore  protected  from  some 
of  the  stability  concerns  which  can  afflict  LQ-type  designs  when  the  dynamic  plant  is  not 
perfectly  known.  Such  stability  concerns  are  usually  very  serious  in  LSS  applications, 
where  one  typically  must  neglect  numerous  known  modes,  in  addition  to  being  forced 
to  neglect  unknown  modes. 

MPC’s  inherent  robustness  enhances  overall  stability  in  the  presence  of  environmen¬ 
tal  uncertainties  or  failed  subsystems.  System  performance,  however,  is  affected  by  the 
accuracy  of  the  on-board  model.  An  on-board  identification  scheme  offers  the  possibility 
to  maintain  good  models,  thereby  positively  impacting  both  performance  and  robust¬ 
ness.  That  is  to  say,  maintenance  of  an  accurate  model  is  maximally  robust,  as  more 
sudden  degradations  can  be  tolerated,  until  system  identification  determines  the  best 
“new”  model. 

At  issue  with  identification  algorithms  is  their  computational  stability,  efficiency  and 
overall  reliability.  The  method  proposed  here,  which  is  based  on  canonical  variate  anal¬ 
ysis,  scores  exceptionally  well  on  these  issues,  has  been  proven  in  related  applications, 
and  is  thus  a  natural  Candida  .e  for  LSS  applications. 

Inclusion  of  an  identification  module  as  a  key  control  subsystem  (Figure  1)  is  recog¬ 
nition  of  the  tight  interdependence  among  modeling,  dynamics  and  control,  as  discussed 
above.  The  estimator  shown  in  Figure  1  is  not  considered  to  be  within  the  scope  of  the 
research  effort,  but  is  recognized  as  being  very  important  in  an  actual  mechanization.  It 
is  assumed  to  be  a  Kalman  filter  or  a  Luenberger-type  reduced  order  observer.  Together 
these  modules  comprise  the  Large  Space  Structure  Identification  and  Control  System 
(LSSICS).  Other  important  modules,  such  as  Failure  Detection  and  Isolation  (FDI),  are 
not  shown  for  purposes  of  clarity. 

LSSICS  has  the  following  attributes  relative  to  crucial  LSS  control  design  issues: 

1.  Identification 

The  CVA  method  of  system  identification  has  b^-en  recently  applied  to  a  variety  of 
systems.  It  has  a  number  of  features  that  are  well  suited  to  real-time  as  well  as  offline 
identification  on  microprocessors  including: 

a)  General  state  space  model  including  inputs,  and  process  and  measure¬ 
ment  noise 

b)  Multi-input  multi-output  systems 

c)  Computation  using  the  singular  value  decomposition 

d)  Numerically  stable  and  accurate  -  never  fails 

e)  Automatic  determination  of  the  best  choice  of  model  state  order 

f)  Finite  amount  of  computation  -  nonrecursive 

g)  No  initialization  required  -  accurate  on  small  samples 

h)  Near  the  maximum  likelihood  lower  bound  in  accuracy 

i)  Discrimination  of  closely  spaced  spectral  peaks 

j)  Simultaneous  identification  of  transfer  function  and  noise  spectrum  due 
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to  disturbance  process 

These  properties  of  the  algorithm  give  it  a  unique  place  among  the  other  currently 
available  algorithms.  It  is  reliable,  accurate,  and  yet  will  handle  the  very  difficult 
problem  of  system  identification  of  multivariable  systems  of  high  order  using  short  data 
lengths. 

2.  Control 

The  proposed  Model  Predictive  Control  method: 

a)  Uses  state  space  form  of  linear  models,  either  as  impulse  or  step 
response  functions 

b)  Is  an  output  predictive  controller.  There  is  thus  no  intrinsic  need 
for  full  state  feedback,  full  order  observers,  etc. 

c)  Like  CVA,  is  structured  to  operate  directly  on  MIMO  systems 

d)  Uses  linear  system  model  in  computing  control  inputs 

e)  Has  demonstrated  robustness  to  model  uncertainties  and/or  changes 

in  the  dynamic  properties  of  the  system  being  controlled . 

Phase  I  of  this  effort  involved  demonstrating  the  feasibility  of  applying  LSSICS  to 
the  LSS  problem.  The  results  obtained  to  date  are  presented  in  Section  4,  and  provide 
further  background  on  the  overall  LSSICS  techniques.  The  main  objective  of  follow- 
on  Phase  II  work  is  to  further  refine  and  time  LSSICS,  integrate  its  identification  and 
control  modules  into  one  efficient  software  entity,  demonstrate  its  capability  on  “real¬ 
istic”  LSS  simulation  models,  compare  its  overall  performance  with  current  alternate 
methods,  and  develop  on-board  software  implementation  requirements. 


2  MATHEMATICAL  FORMULATION 


The  basic  objective  of  the  Phase  I  research  is  to  demonstrr'  V,  'ha  simulation,  the  feasi¬ 
bility  of  applying  MPC  control  and  CVA  id  utification  to  the  design  of  robust  controllers 
for  flexible  space  structures.  Lesser  attention  has  been  given  to  other  key  elements  of  a 
fully  operational  system,  for  example  the  Fault  Detection  and  Isolation  (FDI)  element. 
The  basic  objective  was  achieved  by  defining  an  appropriate  baseline  Large  Space  Struc¬ 
ture  (LSS)  model  and  environment,  and  establishing  test  cases  for  performance  analysis. 
The  MPC  and  CVA  algorithms  were  verified  for  this  new  application,  and  modular  sim¬ 
ulation  tests  were  run  both  for  design  and  performance.  In  this  section,  mathematical 
background  and  formulation  of  the  LSS  model  and  dynamic  process  is  presented,  as 
well  as  salient  features  of  the  CVA  and  MPC  algorithms. 

2.1  LSS  Dynamic  Model 

The  LSS  is  a  distributed  parameter  system  -  consequently,  it  is  infinitely  dimensional 
in  theory,  and  potentially  very  large  in  practice  (analysts  are  working  now  with  models 
containing  10,000  modes!).  However,  implementable  controllers  are  difficult  to  obtain 
for  a  distributed  parameter  system  (Balas,  1982),  a  fact  which  also  applies  to  Model  Pre¬ 
dictive  Control.  For  this  reason,  LSS  dynamics  are  usually  converted  to  a  reduced  order 
model  for  control  synthesis  especially,  but  also  for  verification  purposes  as  a  “truth” 
model.  We  now  summarize  briefly  the  development  of  the  basic  structural  dynamic 
equations. 

A  generic  structure  could  be  described  as  a  continuum  via  the  following  partial 
differential  equation: 


m(x,y,z)uu(x,y,z,t)  +  D0ut(x,y,z,t)  +  A0u{x,y,  z,t)  =  F[x,y,z,t )  (1) 

where  u(x,y,z,t)  is  a  function  of  displacements  of  the  structure  from  its  equilibrium 
position,  which  results  from  the  applied  force  distribution  F(x,y,z,t)  -  also  a  vector  - 
and  from  unmodeled  and  transient  disturbances;  also,  Da  and  Aa  are  linear  operators 
in  x,  y,  z  and  m(z,  y,  z)  is  the  (positive)  mass  density  measure.  The  force  distribution 
vector  and  displacement  vector  may  include  torques  and  rotations,  respectively. 

The  most  popular  method  in  structural  analysis  for  converting  Equation  (1)  to  a 
reduced  order  model  is  the  finite  element  method  (  ref’s  III.3  to  III.5  in  Balas,  1982  ). 
A  set  of  finite  elements  or  patch  functions  5*(x)  is  used  in  approximating  the  solution 
to  Equation  (1)  by 

No 

u{x,t)  =  X  ? k(t)0k(x )  (2) 

Js=l 

for  structural  members  homogenous  in  the  (y,  z)  directions,  where  the  q *  axe  chosen  to 
minimize  the  mean  square  equation  error  when  (2)  is  substituted  into  Equation  (1),  and 
No  is  the  number  of  modes  used  in  the  reduced  order  model.  System  order  N  is  thus 
2 N0.  The  0jk(x)  may  be  linear  combinations  of  piecewise  linear  functions,  cubic  splines 
or  the  actual  LSS  mode  shapes,  when  the  latter  are  known. 
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Defining  displacement  coordinates  q(t)  —  [gl5 q^0}r ,  Equation  (l)  becomes 

Mq  +  Dq  +  Kq  =  B°f  (3) 


The  displacement  equation  (3)  is  readily  converted  into  (linear)  state  space  form  for 
CVA  and  MPC  use;  however,  it  is  common  for  computational  efficiency  and  stability  to 
convert  Equation  (3)  to  modal  form  via  the  orthogonal  transformation 

q  =  Uu 

resulting  in 

u  +  Dit  +  A.u  =  Bf  (4) 

Hence 

x  =  Ax  +  Bf  (5) 

is  the  final  desired  state  space  form,  where  x  =  [urur]T  and 


0  In 

,  B  = 
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and,  considering  also  Equations  (3)  and  (4),  we  have  the  orthogonal  transformations 
(defining  U ) 

UtMU  =  IN ,  UtKU  =  A; 

also 

B  =  UtB°  ,  D  =  UtDU 


and,  finally,  M,D,K,  and  B°  in  Equation  (3)  axe,  respectively,  mass,  damping,  stiffness, 
and  control  influence  matrices.  There  is  a  companion  observation  equation  to  Equation 

(5), 

y  -  Cx  (6) 

As  formulated,  Equations  (5)  and  (6)  are  the  set  of  control  model  equations  used  in 
MPC  design.  Model  error  vectors  representing  disturbances,  unmodeled  modes,  biases, 
etc.,  are  unmodeled  in  the  controller  design  equations  but  generally  used  additively  in 
the  “truth”  model  equivalent  to  Equations  (5)  and  (6).  Structural  analysis  programs 
such  as  NASTRAN  are  generally  used  to  generate  the  reduced  order  model. 

It  is  very  important  that  the  models  derived  by  the  methods  cited  above  be  suffi¬ 
ciently  accurate.  This  fact  is  critical  to  controlling  the  LSS,  and  it  also  motivates  the 
need  for  stable,  reliable  on-board  identification  (see  Section  2.2). 

7 


■  3 .  OOt 


■  6.000  kg 


3  RICIO  BOOY  MOOES 
(ATTITUDE  ANCLES) 

ANO  41  FLEXIBLE  MOOES 
RETAINED  FOR  ACTIVE 
CONTROL  EVALUATION 
(CSCL  MODES  120.  33. 
ANO  41  DISCARDED) 


HEICHT-  21  m 
TOTAL  WEICHT  -  3300  kg 

•  S3  NODES 

•  306  DECREES-OF-FREEDOM 

•  84  NONZERO  MASS  ELEMENTS 

•  127  STIFFNESS  ELEMENTS 


Figure  2:  CSDL  No.  2  Model  (Artist’s  Conception) 


2.1.1  Nonlinear  Models 

Analyses  provided  by  such  programs  as  NASTRAN  are  quite  adequate  for  small  LSS 
motions,  but  become  suspect  as  displacement  amplitudes  increase,  particularly  with 
regard  to  the  size  of  the  rotational  rigid  body  modes.  Large  rigid  body  angles  can 
usually  be  accommodated,  but  if  the  angular  rates  grow  too  much,  certain  nonlinear 
dynamic  effects  have  to  be  modeled,  even  though  structural  deformations  can  still  be 
represented  by  linear  equations.  These  effects  can  be  modeled  by  employing  Grumman ’s 
SATSIM  or  SPACE14  programs  (see  Appendix  A). 

A  version  of  CVA  has  recently  been  postulated  to  identify  certain  general  classes  of 
nonline<uities  (  Larimore,  1987),  but  this  more  complex  version  of  CVA  has  yet  to  be 
fully  validated  in  code.  Similarly,  a  stochastic  version  of  MPC  can  feasibly  be  developed 
to  work  with  nonlinear  models  in  a  manner  compatible  with  CVA  identification;  software 
validation  is  now  underway  on  this  code.  The  focus  of  the  proposed  effort  thus  will 
remain  on  linear  reduced  order  models. 

Although  structural  and  LSS  models  of  varying  degrees  of  dynamic  accuracy  will 
be  used  to  develop  and  adapt  software,  the  primary  baseline  model  at  this  time  is 
expected  to  be  the  CSDL2  structure,  developed  for  the  ACOSS  program  (Strunce,  et 
al.,  1980).  See  Figure  (2).  It  is  expected  that  sufficient  modes  can  be  retained  from 
the  baseline  model  to  display  meaningful  dynamic  “pathologies”  for  proper  analysis  of 
the  identification  and  control  techniques,  and  yet  result  in  a  model  small  enough  for 
efficient  simulation. 


2.1.2  Factors  in  Model  Selection 


A  major  reason  for  pursuing  this  entire  line  of  research  is  our  feeling  that  the  proposed 
identification  and  control  methods  offer  great  potential  for  being  more  cost  effective, 
and  even  outperforming,  currently  available  methods.  This  belief,  as  stated  in  earlier 
sections,  is  based  on  the  excellent  results  obtained  using  CVA  for  identification  in  a 
real-time  wing  flutter  suppression  application,  and  on  the  preliminary  results  we  have 
thus  far  obtained  in  LSS  applications  as  part  of  Phase  I  of  this  effort.  MPC  has  not  been 
developed  to  this  extent,  but  all  of  its  applications  in  a  simulation  environment  point 
to  successful  real-time  LSS  application.  These  applications  include  terrain-following 
(Reid,  et  al.,  1981)  and  battle  damage  reconfiguration  of  aircraft  flight  control  systems 
(Carroll  and  Mahmood,  1986). 

It  is  of  great  use,  then,  to  compare  LSSICS  results  with  those  of  the  other  schemes, 
and  this  comparison  can  be  more  meaningful  if  the  same  baseline  model  is  used.  We 
are  thus  led  to  give  serious  consideration  to  using  a  structure  developed  by  Draper 
Laboratory  for  the  ACOSS  program,  the  CSDL2  model.  This  model  and  its  attributes 
is  detailed  in  Strunce,  et  al.,  (1980).  This  model  represents  a  wide-angle,  three-mirror 
optical  space  system,  about  90  feet  high  and  weighing  10  tons.  Through  the  kindness  of 
Tim  Henderson  and  Dan  Hegg  of  CSDL,  we  have  obtained  a  data  tape  of  NASTRAN 
CSDL2  output,  and  have  developed  code  for  reducing  these  data  to  linear  state  space 
models. 

The  full  model  has  upward  of  50  modes.  Our  Phase  n  work  will  establish  am  ap¬ 
propriate  number  of  modes  for  our  “truth”  model,  and  for  control  design.  The  “truth” 
model  will  perhaps  consist  of  15  to  20  modes,  and  the  control  design  model  at  least  half 
of  that. 

2.1.3  Disturbance  Modeling 

Disturbance  effects  are  not  modeled  directly  in  the  control  synthesis  process,  but  are 
a  part  of  the  “truth”  model.  Additive  white  noise  or  low  order  Markov  processes  will 
be  added  to  the  basic  dynamic  and  measurement  equations.  Random  secular,  or  bias, 
quantities  can  similarly  be  added.  These  would  simulate  certain  environmental  effects 
such  as  solar  pressure,  or  some  types  of  system  failure  -  e.g.,  control  jet  stuck  on. 

We  note  also  that  CVA  identification  also  obtains  models  for  the  above  disturbance 
processes.  This  information  would  be  of  significant  use  to  the  FDI  and  reconfiguration 
functions  (not  addressed  in  this  report)  of  an  LSS  control  system. 

Because  of  the  desire  to  analyze  more  detailed  LSS  models,  and  to  perform  other 
types  of  analyses,  we  realize  the  need  for  greater  computational  resources  than  were 
utilized  in  Phase  I  of  this  project.  We  will  utilize  the  resources  of  our  subcontractor, 
Grumman  Corporate  Research  Center,  to  continue  this  analysis  in  Phase  H. 

2.2  CVA  Identification 

In  this  section,  the  technical  approach  to  LSS  model  identification  is  discussed  in  detail 
including  the  CVA  approach,  selection  of  model  order  and  structure  using  entropy  based 


methods  such  as  the  Akaike  Information  Criterion  (AIC),  and  computational  aspects. 
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2.2.1  Canonical  Variate  Analysis 
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The  generalized  canonical  variate  analysis  approach  has  recently  provided  a  com¬ 
pletely  general  solution  to  the  static  reduced  rank  stochastic  prediction  problem  which 
is  well  defined  statistically  and  computationally  even  when  some  or  all  of  the  various 
covariance  matrices  are  singular  (Larimore,  1986).  All  other  previous  methods  in  the 
statistical  literature  do  not  address  the  general  problem.  For  the  general  time  series 
system  identification  problem,  this  result  guarantees  that  the  solution  to  the  problem  is 
always  well  conditioned  and  produces  a  statistically  meaningful  solution.  In  this  section, 
the  background  for  the  CVA  method  and  details  of  this  new  result  are  given. 

The  analysis  of  canonical  correlations  and  variates  is  a  method  of  mathematical 
statistics  developed  by  Hotelling  (1936;  also  see  Anderson,  1958).  Concepts  of  canoni¬ 
cal  variables  for  representing  random  processes  were  explored  by  Gelfand  and  Yaglom 
(1959),  Yaglom  (1970),  and  Kailath  (1974).  The  initial  application  of  the  canonical 
correlation  analysis  method  to  stochastic  realization  theory  and  system  identification 
was  done  in  the  pioneering  work  of  Akaike  (1974a,  1975,  1976).  This  initial  work  has 
a  number  of  limitations  such  as  no  system  inputs,  no  additive  measurement  noise,  sub¬ 
stantial  computational  burden  involving  numerous  SVD’s,  a  heuristic  set  of  decisions 
for  choosing  a  basis  for  representation  of  the  system,  and  a  number  of  approximations 
including  computation  of  the  AIC  criterion  for  decision  on  model  order. 

Some  important  generalizations  and  improvements  in  Akaike ’s  canonical  correlation 
method  have  recently  been  made  by  Larimore  (1983b).  These  include  generalization  to 
systems  with  additive  measurement  noise  and  with  inputs  including  feedback  controls. 
A  major  departure  of  the  approach  from  previous  work  is  the  use  of  a  single  SVD  to 
optimally  choose  k  linear  combinations  of  the  past  for  prediction  of  the  future.  The  very 
natural  measure  of  quadratically  weighted  prediction  errors  at  possibly  all  future  time 
steps  is  used.  Formulated  as  such  a  prediction  problem,  it  is  shown  how  a  generalized 
canonical  variate  analysis  gives  the  solution  explicitly.  The  interpretation  of  canonical 
variates  as  optimal  predictors  is  central  in  motivating  interest  in  such  a  problem  formu¬ 
lation  and  is  scarcely  found  in  the  statistical  literature  (Larimore,  1986).  The  optimal 
k-order  predictors  are  not  in  general  recursively  computable,  but  the  optimal  state- 
space  structure  for  approximating  them  is  expressed  simply  in  terms  of  the  canonical 
variate  analysis.  The  problem  of  finding  an  optimal  Hankel  norm  reduced  order  model 
(Adamjan  et  al,  1971;  Kung  and  Lin,  1981)  is  related  to  the  canonical  variate  approach 
(Camuto  and  Menga,  1982;  Larimore,  1983b).  The  balanced  realization  method  is  a 
particular  case  of  the  generalized  canonical  variate  analysis  (Desai  and  Pal,  1984).  To 
more  concisely  discuss  the  related  research,  the  problems  of  identification,  reduced  order 
modeling  and  filtering  can  be  described  as  follows  (Larimore,  1983b). 

Consider  the  problem  of  choosing  an  optimal  system  or  model  of  specified  order 
for  use  in  predicting  the  future  evolution  of  the  process.  Consider  the  past  vector  pt 
consisting  of  past  outputs  yt  and  inputs  ut  before  time  t  and  the  future  vector  ft  of 


outputs  at  time  t  or  later  so 


pj  =  (yt-i,yi-j,  fj  =  (yt,yt+i,---)T  (7) 

We  assume  that  the  processes  yt  and  ut  are  jointly  stationary  and  denote  the  covariance 
matrices  among  /  and  p  as  E//,EPP,  and  E/p. 

The  major  interest  is  in  determining  a  specified  number  k  of  linear  combinations 
of  the  past  pt  which  allow  optimal  prediction  of  the  future  ft.  The  set  of  k  linear 
combinations  of  the  past  pt  are  denoted  as  a  k  x  1  vector  mt  and  are  considered  as 
Ar-order  memory  of  the  past.  The  optimal  linear  prediction  ft  of  the  future  /*,  which  is 
a  function  of  a  reduced  order  memory  mt,  is  measured  in  terms  of  the  prediction  error 

£{ll  /.  -  h  III,}  =  E{(f,  -  f,)TL'(h  -  /,)>  (8) 

where  E  is  the  expectation  operation  and  A  is  an  arbitrary  positive  semidefinite  sym¬ 
metric  matrix  so  that  A^  is  an  arbitrary  quadratic  weighting  that  is  possibly  singular. 
The  optimal  prediction  problem  is  to  determine  an  optimal  fc-order  memory 


mt  =  Jkpt  (9) 

by  choosing  the  k  rows  of  Jk  such  that  the  optimal  linear  predictor  based  on  mt 

minimizes  the  prediction  error  (8). 

As  derived  in  Larimore  (1986),  the  solution  to  this  problem  in  the  completely  general 
case  where  the  matrices  E //,  Ep,,,  and  A  may  be  singular  is  given  by  the  generalized 
singular  value  decomposition  as  stated  in  the  following  theorem. 

Theorem  1.  Consider  the  problem  of  choosing  k  linear  combinations  mt  =  Jkpt 
of  pt  for  predicting  ft,  such  that  (8)  is  minimized  where  Epp,  and  A  are  possibly  sin¬ 
gular  positive  semidefinite  symmetric  matrices  with  ranks  m  and  n  respectively.  Then 
the  existence  and  uniqueness  of  solutions  are  completely  characterized  by  the  (EPP,A)- 
generalized  singular  value  decomposition  which  guarantees  the  existence  of  matrices  J, 
L,  and  generalized  singular  values  71,  •  •  •  ,7r  such  that 

JYjppJ  =  Im,  LAX  =  /n,  JEP/X  =  Diag(~fi  ^  ^  7r  ^  0,  ■  •  • ,  0)  (10) 

The  solution  is  given  by  choosing  the  rows  of  Jk  as  the  first  k  rows  of  J  if  the  A:-th 
singular  value  satisfies  7*  >  7*+1.  If  there  are  r  repeated  singular  values  equal  to  7*, 
then  there  is  an  arbitrary  selection  from  among  the  corresponding  singular  vectors,  i.e. 
rows  of  J.  The  minimum  value  is 

h  ~  h  &>  =  *ri,E"  -*-■—&  <»> 

This  result  not  only  gives  a  complete  characterization  of  the  solutions  in  selecting 
optimal  predictors  mk  from  the  past  pt  for  prediction  of  the  future  /t,  but  the  reduction 
in  prediction  error  for  all  possible  selections  of  order  k  is  given  simply  in  terms  of  the 
generalized  singular  values.  This  is  of  great  importance  since  it  avoids  having  to  do  a 
considerable  amount  of  computation  to  determine  what  selection  of  order  is  appropriate 
in  a  given  problem. 


Different  selections  of  the  weighting  matrix  A  can  be  used  for  different  purposes. 
A  number  of  classical  reduced  rank  statistical  analysis  problems  of  static  variables,  i.e. 
with  independence  from  ‘time’  to  ‘time’,  can  be  formulated  and  solved  by  the  generalized 
CVA  of  Theorem  1.  In  the  classical  canonical  correlation  analysis  problem,  A  =  E//.  In 
the  principal  components  analysis  problem,  the  ‘past’  and  ‘future’  are  the  same  space 
and  in  addition  A  =  I.  A  generalization  of  this  with  the  ‘past’  and  ‘future’  different  is 
the  principal  component  analysis  of  instrumental  variables  where  A  =  I  (Rao,  1965). 
The  only  consideration  of  the  case  of  singular  covariance  matrices  for  the  canonical 
correlation  analysis  problem  is  by  Khatri  (1976).  The  solution  is  considerably  more 
complicated  and  is  not  related  to  a  computational  procedure  such  as  the  SVD.  Also  it 
does  not  address  the  general  CVA  problem  with  A  an  arbitrary  positive  semidefinite 
matrix.  For  system  identification,  the  use  of  the  weighting  matrix  A  =  E {j  results  in 
a  near  maximum  likelihood  system  identification  procedure  (Larimore,  Mahmood  and 
Mehra,  1984). 

2.2.2  Selection  of  State  Space  Order  and  Model  Structure 

The  generalized  CVA  method  allows  the  determination  of  the  fit  of  the  various  state 
space  models  and  the  selection  of  the  best  model  state  order  before  computation  of  the 
state  space  models.  This  “built-in”  feature  of  the  CVA  method  has  great  significance 
with  regard  to  the  LS8  problem,  for  which  it  may  be  quite  risky  merely  to  truncate 
high  frequency  modes.  The  model  order  selection  problem  has  several  aspects: 

(1)  the  selection  of  the  best  model  order  in  some  statistical  sense  based  upon 
the  observed  data,  and 

(2)  model  order  reduction  based  upon  some  criterion  such  as  control  performance 
error,  and  not  merely  magnitude  of  modal  frequency. 

Both  of  these  issues  are  addressed  by  approaching  the  problem  using  the  generalized 
canonical  variate  analysis.  In  the  remainder  of  this  section,  the  system  identification 
problem  is  discussed  including  the  statistical  order  determination  and  order  reduction 
problems. 

Consider  the  general  case  of  the  reduced  order  filtering  and  modeling  problem:  given 
the  past  of  the  related  random  processes  ut  and  yt,  we  wish  to  model  and  predict  the 
future  of  yt  by  a  fc-order  state  xt  and  state-space  structure  of  the  form 

X{+i  =  F  Xt  +  Gut  +  Wt  (12) 

yt  =  Cxt  +  Aut  +  Bwt  +  vt  (13) 

where  xt  is  the  state  and  wt  and  vt  are  white  noise  processes  that  are  independent  with 
covariance  matrices  Q  and  R  respectively.  Equations  (12)  and  (13)  may  be  seen  as 
discretizations  of  Equations  (5)  and  (6),  with  allowance  here  for  feedforward  dynamics 
and  random  processes.  The  matrices  A  and  B  have  a  different  context  here  than  in 
(5).  The  white  noise  processes  model  the  covariance  structure  of  the  error  in  predicting 
yt  and  x*+i  from  ut  and  xt.  A  special  case  of  the  reduced-order  filtering  problem  is 

the  transfer  function  approximation  problem  where  ut  and  yt  are  the  input  and  output 

processes  and  an  approximate  state-space  model  is  desired. 
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In  the  computational  problem  given  finite  data,  the  past  and  future  of  the  process 
are  taken  to  be  finite  of  length  d  lags  so 

p!  =  (vf-1,  •  •  • ,  yf-d,  «r_i,  •  •  • ,  «r_  d)T,  fj  =  (yf ,  •  •  • ,  yJ+d-\)T  (14) 

Akaike  (1976)  proposed  choosing  the  number  d  of  lags  by  least  squares  autoregressive 
modeling  using  recursive  least  squares  algorithms  and  choosing  the  number  of  lags 
as  that  minimizing  the  AIC  criterion  discussed  below.  This  insures  that  a  sufficient 
number  of  lags  are  used  to  capture  all  of  the  statistically  significant  behavior  in  the  data. 
This  procedure  is  easily  generalized  to  include  the  case  with  inputs  it.  In  the  model 
identification  problem,  by  using  the  weighting  matrix  A  =  2//  the  identified  system  is 
close  to  the  maximum  likelihood  estimation  solution  (Lari more,  Mahmood,  and  Mehra, 
1984).  The  generalized  SVD  of  Theorem  1  above  determines  a  transformation  J  of  the 
past  that  puts  the  state  in  a  canonical  form  so  that  the  memory  mf  =  Jkpt  contains 
the  states  ordered  in  terms  of  their  importance  in  modeling  the  process.  The  optimal 
memory  for  a  given  order  k  then  corresponds  to  selection  of  the  first  k  states. 

For  the  determination  of  model  state  order,  recent  developments  in  the  selection 
of  model  order  and  structure  based  upon  entropy  or  information  will  be  used.  Such 
methods  were  originally  developed  by  Akaike  (1973)  and  involve  the  use  of  the  Akaike 
Information  Criterion  (AIC)  for  deciding  the  appropriate  order  of  a  statistical  model. 
The  AIC  for  each  state  order  k  is  defined  by 

AIC(k)  =  -2  log  p(Y",  UN ;  0k)  +  2 Mk  (15) 

where  p  is  the  likelihood  function,  based  on  the  observations  (YN,  UN)  at  N  time  points, 
with  the  maximum  likelihood  parameter  estimates  9k  using  a  fc-order  model  with  Mk 
parameters.  The  model  order  k  is  chosen  with  the  minimum  value  of  AIC(fc).  A 
predictive  inference  justification  of  the  use  of  an  information  or  entropy  based  criterion 
such  as  AIC  is  given  in  Larimore  (1983a)  based  upon  the  fundamental  principles  of 
sufficiency  and  repeated  sampling.  The  number  of  parameters  Mk  in  the  state  space 
model  (12)  and  (13)  is  determined  by  the  general  state  space  canonical  form  as  in  Candy 
et  al  (1979)  as 

Mk  =  2  kn  +  km  +  nm  +  n(n  +  l)/2  (16) 

where  n  and  m  are  the  dimensions  of  the  output  and  input  vectors  yt  and  ut  respectively. 
This  is  far  less  than  the  number  of  elements  in  the  various  state  space  matrices. 

In  the  paper  of  Akaike  (1976),  an  approximate  procedure  for  evaluating  the  AIC 
using  the  sample  canonical  correlations  is  given.  This  procedure  has  been  found  to  be 
highly  approximate  because  the  canonical  correlation  theory  upon  which  it  is  based  as¬ 
sumes  independence  of  variables  which  is  violated  for  a  correlated  time  series.  The  exact 
AIC  can  be  evaluated  directly  from  the  generalized  CVA  described  above,  although  it 
requires  substantial  computation.  Recently,  efficient  algorithms  for  accurate  evaluation 
of  the  AIC  from  the  generalized  CVA  have  been  developed  at  BTS. 

In  choosing  model  order,  the  risks  of  introducing  bias  into  the  model  in  choosing 
too  low  an  order  must  be  weighed  against  introducing  additional  variability  in  choosing 
too  high  an  order  (Larimore  and  Mehra,  1985).  Having  selected  the  optimal  statistical 
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order  for  the  model  of  the  observed  data,  additional  order  reduction  may  be  appropriate 
in  determining  a  filter  or  controller  in  various  applications.  As  discussed  in  Larimore 
(1983b),  the  general  weighting  criterion  used  in  the  generalized  CVA  gives  a  natural 
procedure  for  such  additional  order  reduction.  This  is  very  important  in  control  design 
for  adaptive  control.  A  general  weighting  can  be  used  which  reflects  the  important 
frequencies  which  are  necessary  to  control. 

Once  the  optimal  fc-order  memory  mt  is  determined,  state-space  equations  of  the 
form  (12)  and  (13)  for  approximating  the  process  evolution  are  easily  computed  by  a 
simple  multiple  regression  procedure  (Larimore,  1983b).  Since  the  CVA  system  identi¬ 
fication  procedure  involves  the  state  space  model  form,  it  has  the  major  advantage  that 
the  model  is  globally  identifiable  so  that  the  method  is  statistically  well-conditioned  in 
contrast  to  ARMA  modeling  methods  (Gevers  and  Wertz,  1982).  Furthermore,  since 
the  computations  are  primarily  a  SVD,  the  computations  are  numerically  stable  and 
accurate  with  an  upper  bound  on  the  required  computations.  Thus  the  method  is  com¬ 
pletely  reliable,  and  has  been  demonstrated  as  such  in  the  adaptive  flutter  suppression 
system  in  wind  tunnel  tests  involving  reidentification  of  the  system  dynamics  tens  of 
thousands  of  times.  From  the  theory  of  the  CVA  method  (Larimore,  Mahmood  and 
Mehra,  1984),  it  can  be  shown  that  there  are  no  difficulties  such  as  biased  estimates 
caused  by  the  presence  of  a  correlated  feedback  signal. 

2.2.3  Computational  Aspects 

For  complex  space  structures,  the  dimension  of  the  matrices  can  easily  be  in  the  hun¬ 
dreds  or  larger.  The  computation  of  the  SVD  for  such  large  matrices  can  require  con¬ 
siderable  time  on  conventional  serial  computers.  The  availability  of  multiprocessors 
such  as  the  CRAY  or  the  INMOS  Transputer  makes  possible  considerable  reduction  in 
the  required  computation  time.  We  propose  here  to  work  with  models  sized  for  the 
Micro- VAX. 

A  very  efficient  algorithm  for  computing  the  usual  SVD  has  been  recently  devised 
for  highly  parallel  systolic  arrays  by  Brent  and  Luk  (1985).  This  algorithm  was  recently 
generalized  (Luk,  1985)  in  one  particular  way  to  compute  the  B-SVD  (Vein  Loan,  1976) 
which  is  different  from  the  generalized  SVD  required  for  the  CVA  method.  An  n  x  n 
square  systolic  array  of  processors  requires  communication  only  between  the  nearest 
neighbor  processors  in  synchrony  with  the  computational  cycle  of  all  the  processors. 
The  computation  of  the  SVD  of  a  nxn  matrix  using  a  nxn  array  of  processors  requires 
only  order  n  processor  cycles  as  compared  to  order  n  cubed  for  a  serial  computer  with 
a  single  processor.  Such  parallel  processors  and  algorithms  could  make  routine  the 
analysis  of  very  large  sets  of  variables  such  as  arise  naturally  in  the  canonical  variate 
analysis  of  large  order  or  nonlinear  systems. 

Modern  computer  algorithms  (Golub,  1969)  for  canonical  correlation  analysis  use  a 
standard  singular  value  decomposition  to  compute  the  generalized  singular  value  decom¬ 
position  (see  Equation  (10))  with  A  =  E//  by  first  finding  square  root  factors  of  Ep,,  and 
A,  and  then  doing  a  standard  singular  value  decomposition  on  A  =  Epp/2Ep/(A-1/2)  = 
QS RT  where  QQT  —  I  =  RRT  and  S  is  diagonal.  Then  the  generalized  singular  value 
decomposition  is  given  by  J  =  QTEpp/2,  L  =  RTA-1^  and  D  =  S.  Thus  the  joint 


orthonormalization  of  pt  and  ft  in  the  norms  £pp  and  A  to  give  the  canonical  covariance 
structure  D  is  very  naturally  viewed  as  a  generalized  singular  value  decomposition  both 
in  terms  of  the  simple  reduction  discussed  in  Theorem  1  as  well  as  the  actual  computa¬ 
tional  algorithms.  This  can  be  determined  computationally  using  a  standard  singular 
value  decomposition  which  is  numerically  very  accurate  and  stable  as  compared  with 
the  earlier  eigenvalue  computational  methods  (Bjorck  and  Golub,  1973).  An  open  topic 
is  the  investigation  of  numerical  methods  that  directly  compute  the  particular  gener¬ 
alized  SVD  rather  than  transforming  the  problem  to  the  standard  SVD.  Such  a  direct 
approach  may  have  better  overall  numerical  accuracy. 

A  second  problem  is  specified  in  terms  of  the  observed  data  given  as  N  observations 
(pi, . . .  ,p//)  =  B  and  (/i, . . . ,  fs)  =  C  on  p  and  /  respectively.  The  usual  sample  covari¬ 
ances  are  computed  as  £pp  =  BBT,  2P/  =  BCT  and  2//  =  CCT  which  mathematically 
are  used  in  the  generalized  singular  value  decomposition.  Numerically,  however,  the 
formation  of  these  products  defining  the  sample  covariances  results  in  a  halving  of  the 
numerical  precision  of  the  computation.  In  the  case  of  given  data,  Bjorck  and  Golub 
(1973)  give  computational  procedures  that  avoid  these  squaring  operations  and  operate 
directly  on  the  observed  data. 


2.3  Model  Predictive  Control 


A  critical  need  for  a  robust,  adaptive  controller  is  a  reliable  system  model  of  appropriate 
order.  We  assume  that  an  appropriate  identification  method  exists  to  supply  these 
models  to  the  controller.  Reduced  order  model  predictive  control  gains  can  be  designed 
adaptively  in  real  time  as  the  models  are  updated.  Control  synthesis  can  thus  take  place 
by  rote  (pre-stored  algorithms),  without  intensive  engineering  analysis,  as  the  system 
model  changes.  A  form  of  model  predictive  control  which  involves  the  use  of  the  Singular 
Value  Decomposition  (SVD)  has  been  outlined  by  Reid  et  al.  (1981),  and  applied  to 
reconfigurable  flight  control  systems  by  Carroll  et  al.,  (1986).  Because  SVD  and  MPC 
are  similar  computationally,  there  is  an  opportunity  for  a  more  efficient  algorithmic 
structure,  especially  if  systolic  array  processors  are  utilized.  This  particular  type  of 
controller  is  thus  very  compatible  with  computationally  stable  identification  methods 
(Section  2.2)  and  is  now  briefly  described. 

The  Model  Predictive  Control  (MPC)  technique  has  been  developed  as  the  next 
generation  of  a  widely  used  control  technique  known  as  Model  Algorithmic  Control 
(MAC).  There  has  lately  been  a  growing  interest  in  the  class  of  MAC-type  controllers, 
known  also  as  “predictive  controllers”,  and  MAC  is  probably  the  earliest  one  in  this 
class  reported  in  the  literature,  as  noted  above. 

MPC  has  its  origins  in  adaptive  controllers  developed  for  industrial  processes.  There 
is  an  input/output  structure  to  this  type  of  controller  which  is  amenable  to  adaptation 
to  system  parameter  changes,  without  the  need  to  understand  explicitly  the  dynamics 
behind  the  particular  change.  This  input/output  structure  eliminates  the  formal  need 
for  full  state  feedback  required  by  most  other  multivariable  control  synthesis  methods, 
and  allows  for  greater  structural  compatibility  with  robust  identification  methods.  An¬ 
other  serious  limitation  of  the  earlier  controllers  is  that  they  typically  use  am  impulse 
response  model  of  the  plant  to  be  controlled  (Richalet,  et  al.  1978)  and  therefore  cannot 
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be  used  for  an  open- loop  unstable  plant. 

The  MAC  and  MPC  techniques  are  based  on  the  same  principles.  These  techniques 
are  fundamentally  and  philosophically  different  from  “feedback”  controllers  such  as  the 
LQ  regulator  and  its  variants,  in  which  there  is  an  explicit  notion  of  “feedback”  of  the 
current  state  to  derive  a  closed-loop  control  law.  Instead,  closed  loop  control  is  achieved 
by  accomplishing,  at  every  cycle  of  the  digital  control  loop,  the  following  functions:  (a) 
predict  the  future  zero  input  response  of  the  plant  using  an  observer  or  a  model  of  the 
plant  over  a  selected  “horizon  of  prediction”.  The  duration  of  this  prediction  interval  is 
a  key  stability  and  robustness  parameter.  The  prediction  must  be  closed  loop,  i.e.,  based 
on  available  measurements  of  the  actual  output  response.  A  Kalman  filter  or  observer 
can  provide  an  estimate  of  the  current  state,  based  on  measured  outputs  and  the  system 
mathematical  model,  (b)  Compute  a  desired  future  reference  trajectory  along  which  the 
system  is  desired  to  be  moved  over  the  “horizon  of  prediction”  to  a  set  point  (which  can 
be  zero  for  a  regulator  problem,  or  a  function  of  time  for  tracking,  e.g.,  terrain  following 
or  on-orbit  target  tracking),  (c)  Find  a  control  sequence  to  minimize  the  Euclidean 
distance  between  the  predicted  output  and  the  desired  trajectory.  Reid,  et  al.,  (1981) 
have  developed  a  way  to  generate  the  future  control  sequence  by  formulating  the  error 
criteria  as  a  linear  least  squares  problem.  This  is  done  by  adding  “output  smoothing 
points”  -  i.e.,  by  sampling  outputs  at  a  faster  rate  than  that  at  which  the  control  is 
updated.  There  is  thus  an  implicit  control  of  the  system  state  even  though  there  is 
no  explicit  state  control.  The  result  is  that  the  output  trajectory  remains  close  to  the 
reference  trajectory  even  between  control  update  times;  this  has  an  effect  of  reducing 
system  internal  energy,  which  for  the  LSS  is  directly  related  to  large  excursions  of  the 
state  variables  (displacement  or  modal  coordinates)  from  their  reference  level  (Colson, 
1978). 

The  difference  between  MAC  and  MPC  lies  in  the  fact  that  MAC  uses  an  impulse 
response  model  of  the  plant  to  predict  its  future  output  and  hence  cannot  be  used  to 
control  an  open-loop  unstable  plant.  To  overcome  this  problem,  a  special  formulation 
which  uses  a  state-space  model  of  the  plant,  MPC,  has  been  developed.  However, 
the  main  advantage  of  MPC/MAC  is  that  there  is  a  clear  and  transparent  relationship 
between  system  performance  and  various  parameters  embedded  in  the  design  procedure. 
This  feature  is  particularly  useful  in  LSS  applications,  where  effective  reduced  order 
modeling  is  vital.  SVD  identification  clearly  enhances  this  feature  as  well.  It  should 
finally  be  pointed  out  that  the  structure  of  MPC  allows  for  decentralized,  or  distributed 
control  designs.  The  peculiarities  of  LSS  dynamics  often  dictate  such  a  design  (Kosut 
and  Lyons,  1984). 

In  this  section  we  develop  an  exact  model  of  the  algorithm  for  control  computa¬ 
tion,  find  an  explicit  form  of  the  control  law,  and  place  MPC  in  the  standard  control 
design  framework  so  that  stability,  gain  and  phase  margin,  robustness,  etc.,  can  be  as¬ 
certained  a  prion’.  The  designer  can  then  iterate  on  various  design  parameters  to  meet 
the  specifications  of  the  particular  LSS  system. 

There  are  many  variations  of  predictive  controllers  such  as  DMC  (Cutler  and  Rrt- 
maker,  1979/80),  IMC  (Garcia  and  Moyari,  1982),  etc.,  but  here  we  closely  follow  the 
developments  in  Reid  et  al.  (1981)  and  Carroll  and  Mahmood,  (1986). 


2.3.1  Basic  Elements  of  MPC 


There  are  five  elements  in  MPC: 

(i)  A  dynamic  plant  to  be  controlled: 

x(fc  +  l)  =  Fx(k)  +  Gu(k),  x(0)  =  xQ  (17) 

y{k)  =  Cx(k)  (18) 

We  assume  that  x(k),u(k)  and  y(k )  are  of  dimension  n,  m  and  p,  respectively.  Also  F 

and  G  axe  the  discrete  system  and  control  matrices,  respectively.  A  feedforward  term  is 
straightforward  to  add  in  Equation  (18),  but  is  avoided  here  for  clarity.  It  is  not  vital 
that  the  plant  model  be  linear. 

(ii)  An  internal  model  of  the  plant  having  the  same  input-output  dimension  m  x  p 
as  that  of  the  best  representation  of  the  actual  plant. 

x(k  +  1)  =  Fx(k)  +  Gu(k),  x(0)  =  x0  (19) 

y(k)  =  Cx(k )  (20) 

(iii)  A  p-dimensional  reference  trajectory  yr(k),  along  which  the  system  is  desired  to 
be  guided  to  a  set  point.  There  are  many  variations  of  generating  this  desired  trajectory. 
In  our  approach,  yr(fc)  is  a  preferably  smooth  curve  initialized  on  the  current  output  of 
the  actual  plant  y(k )  that  leads  y{k)  to  a  possibly  time  varying  p-dimensional  set  point 
s(k).  In  this  analysis  yr{k)  evolves  on: 

yr(fc  +  l)  =  A.ayr(k)  +  (I  —  Aa)s(k),  yr(k)  =  y(k);  (21) 

|A,(Aa)|  <  1,  for  all  i, 

where  At  (X)  is  the  \th  eigenvalue  of  X. 

Usually  Aa  is  a  diagonal  matrix,  i.e.,  Aa  =  diag(oti)  and  0  <  a,  <  1.  The  reference 
trajectory  in  Equation  (21)  has  first  order  dynamics  but  may  not;  splines  or  some  other 
mechanism  of  generating  the  reference  trajectory  can  be  used.  However,  it  is  easy  to 
see  that  the  solution  of  (21)  is,  if  s(fc)  =  s, 

yr(fc  +  lV)=A"y(fc)  +  (/  — A")s  (22) 

The  a,  thus  are  like  discrete  first  order  time  constants  which  can  be  adjusted  for  a 
certain  accuracy  in  yr  after  N  steps.  Any  method  chosen  should  be  simple,  as  the 
calculation  is  done  every  control  update. 

(iv)  A  closed-loop  prediction  scheme  for  predicting  the  future  output  of  the  plant 
according  to 


y p{k  +  1)  =  y{k  +  1)  +  y{k)  -  y{k)  (23) 

The  yp(k  +  1)  may  be  considered  to  be  elements  of  Y/ . 

Again,  we  have  used  this  scheme  for  the  sake  of  simplicity  but  the  user  can  supply 
his  own  routine  for  generating  the  predictions.  An  acceptable  implementation  is  to  use 
the  zero  input  response  to  the  state  at  k,x(k)  :  CF'x(k),  i  =  1  ,...,iV,  for  each  of  the 
N  outputs  over  the  horizon  of  prediction. 
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(v)  A  quadratic  cost  functional  J  based  on  the  error  between  yp(k)  and  yr(k )  over  a 
finite  horizon  N: 

N 

J  =  tr  Y]\Q(k  +  i)e(k  +  J )eT(k  +  f)  +  R(k  +  i  -  l)it(A:  +  i  —  l)ur(fc  +  t  —  1)]  (24) 

«=i 

where  Q(-)  and  R(- )  are  positive  semi-definite  possibly  time-varying  weights  and  typi¬ 
cally  e(k  +  i)  is  the  error  to  be  minimized  in  the  overdetermined  problem  (see  below). 
In  most  applications  R(-)  is  set  to  zero.  This  is  acceptable  since  physical  limits  on  the 
control  element  excursions  are  easily  incorporated  into  the  logic  for  deriving  the  control 
sequence. 

Given  (i)  -  (v),  MPC  finds  an  optimal  control  sequence  u'(k  +  i  —  1),  i  =  1, 2, ...,  N 
by  minimizing  J  over  the  admissible  input  sequence  u(k  +  i  —  1)  e  Cl,  t  =  1,2, 

Once  this  sequence  is  computed,  the  first  element,  i.e.,  ti'(k),  is  applied  to  the  actual 
plant  and  the  entire  process  is  repeated  all  over  again.  Note  that  the  input/output 
nature  of  MPC,  and  specifically  the  direct  incorporation  of  a  predicted  desired  output 
in  the  calculation  of  the  optimal  control  sequence,  mean  that  merely  a  simple  structuring 
of  the  output  setpoints  s(fc)  will  cause  MPC  to  behave  as  a  regulator  (shape  control)  or 
a  tracker  (vibration  control). 

Usually  a  control  law  is  designed  on  the  assumption  that  the  actual  plant  is  the  same 
as  the  nominal  model  although  the  latter  is  almost  invariably  different  from  the  former. 
The  robustness  of  the  designed  control  law  is  evaluated  on  the  nominal  model  which 
specifies  the  stability  neighborhood  around  it.  If  the  actual  plant  lies  in  this  region  the 
nominal  control  law  will  guarantee  the  closed-loop  stability  for  the  actual  plant. 

As  Equation  (24)  implies,  MPC  control  can  be  formulated  as  a  weighted  least  squares 
problem,  which  develops  as  follows2: 

Future  output  y(k)  is  related  to  present  state  i(0)  and  future  inputs  (u(0), u(l), ..., 
u(k  —  1)}  via  the  discrete  equation 

y{k)  =  yti(k)  +  y„(k)  (25) 

(Note  that  current  time  is  referenced  to  k  =  0.)  In  Equation  (25),  is  the  zero  input 
response, 

yn(k)  =  CFkx(  0), 

and  y„  is  the  zero  state  response,  or  (see,  e.g.,  Reid  1983), 

y*>ik)  =  E  Mfc  “  *>(*’)  (26) 

«=o 

This  formulation  assumes  that  the  control  input  is  piecewise  constant  over  the  control 
sample  interval,  Te ,  and  indexed  so  that 

u(f)  =  u(k),  t  e{kTe,(k +  1)TC) 

2Thia  development  is  based  on  Reid,  et  al.,  (1981),  but  is  expanded  to  handle  MIMO  problems. 
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We  now  establish  other  preliminaries:  let  NSM  be  the  number  of  “output  smoothing” 
points  per  control  update.  We  assume  throughout  this  rep  .rt  that  the  system  dis¬ 
cretization  and  output  sampling  periods,  Td,  are  the  same.  This  assumption  is  one  of 
convenience  which  may  be  relaxed  at  a  future  time  without  much  difficulty,  if  warranted. 
Thus,  Tc  —  NSM  •  Td-  If  L  is  the  number  of  discrete  controls  to  be  predicted  into  the 
future,  then  the  total  number  of  samples  in  the  prediction  horizon  is  NSM  •  L.  It  turns 
out  that  NSM  and  L,  in  addition  to  the  cost  weighting  matrices  Q  and  R,  are  key  MPC 
design  parameters. 

Equation  (26)  represents  a  discrete  convolution  process,  with  hd(k)  being  the  discrete 
impulse  response  function 

hd{k)  =  C  Fk~lG  (27) 

The  MPC  problem  then  involves  solving  Equation  (24)  for  the  sequence  (u(fc)},  taken 
over  the  horizon  of  prediction,  with  y(k)  in  (24)  replaced  by  the  predicted  desired 
output,  t id{k)-  This  is  done  by  considering  the  NSM  •  L  outputs  in  the  horizon  of 
prediction  simultaneously,  substituting  Equation  (26)  into  (25),  and  formulating  as  a 
vector-matrix  equation 

Yzi  +  HU  =  Yd  (28) 

We  note  here  that  the  Y-terms  in  Equation  (28)  consist  of  a  string  of  NSM  •  L  p- 
dimensional  output  vectors: 


Mi) 

(i) 

y<<(2) 

..  y»(2) 

Yd  = 

,  Yzi  =  . 

.  y„(NSM  •  L)  . 

[  y„(NSM  •  L)  . 

The  yd{k)  are  generated  according  to  designer  option,  performance  requirements,  etc. 
(see,  for  example,  Equation  (21)).  This  is  usually  done  according  to  whether  the  control 
objective  is  to  achieve  some  constant  value  (“setpoint”)  s,  the  regulator  case,  or  to 
follow  a  varying  reference  path,  the  tracking  case.  The  LSS  problem  has  examples  of 
both  general  cases:  regulation  of  shape  (eg.,  antenna)  and  step  changes  in  attitude 
are  examples  of  a  regulator  mode,  while  target  tracking  is  obviously  a  tracker  mode 
example. 

The  explicit  formulation  of  HU  depends  on  such  issues  as  the  desirability  of  using  im¬ 
pulse  response  or  step  response  functions,  etc.  Since  these  variants  are  mathematically 
equivalent,  the  practical  issues  are  then  computational  accuracy,  availability,  storage 
and  run-time  efficiency  -  tradeoffs  which  we  propose  to  analyze  in  this  project.  A  fa¬ 
vored  mechanization  is  to  eliminate  some  summations  by  using  step  response  functions 
instead  of  impulse  response  functions: 

r(k)  =  Y  MO  =  r(k  “  0  +  Mfc) 


where  r(k )  is  the  step  response  at  sample  point  k  to  a  unit  step  input  applied  at  reference 
sample  point  k  =  0.  Rewriting  Equation  (26)  thus  means  that 


U  4 


u(0) 

r  ^ 

0 

••  O' 

Au(l) 

h2 

hi  • 

••  0 

: 

,  H  = 

: 

: 

A  u(L  —  1) 

.  hL 

hL-i  • 

■  ^  . 

(29) 


where  the  Au(i)  are  changes  to  the  previous  applied  control,  i.e. , 

u(k)  =  u(0)  +  Au(l)  + ...  +  A  u(k)  =  u(k  —  1)  +  A  u(k) 

U  is  an  ml  x  1  vector  array.  H  is  very  similar  in  form  (and  information  content)  to 
the  Hankel  matrix.  The  h*  in  H  are  submatrices  consisting  themselves  of  NSM  transfer 
function  matrices  each  of  size  p  x  m.  Dimension  (ht)  =  p  •  NSM  x  m  in  the  following: 


hit  = 


'  r(NSM(fc  -  1)  +  1)  ‘ 

r(NSM  ■  k) 

(30) 


where  the  hd  are  defined  in  Equation  (27).  Equation  (30)  is  a  corrected  version  of 
Equation  (15)  in  Reid,  et  al.,  (1981),  and  it  is  extended  to  the  MIMO  case,  as  are  the 
other  equations  in  this  development. 

Having  formulated  the  MPC  problem,  the  solution  is  to  solve  Equation  (28)  for  U. 
Reid’s  contribution  was  to  simultaneously  control  internal  plant  energy  and  convert 
Equation  (28)  into  an  overdetermined  linear  equation  problem  by  adding  the  NSM 
output  smoothing  points  per  control  cycle  (Reid,  et  al.,  1981).  The  usual  least  squares 
“approximate”  solution  methods  then  apply.  The  weighted  criterion  is  needed  to  keep 
(28)  as  valid  as  possible,  and  also  (an  extension  of  Reid,  developed  by  him  and  Carroll 
and  Mahmood  (1986))  to  penalize  large  control  excursions.  The  cost  criteri  (24)  then 
becomes 


J  =  [HU  -  [Yd-Yzi)\T  Q[HU  -  [Yd-Yzi)\  +  UT RU  (31) 

which  has  for  a  solution 

U"  =  [HtQH  +  R)~1HrQ{Yd-Yzi)  (32) 

Note:  (l)  Equation  (32)  solves  simultaneously  for  the  L  future  controls  (cf.,  Equation 
(29)),  but  only  the  first  element,  u(0),  is  actually  applied  to  the  system;  in  a  closed  loop 
operation  the  problem  will  be  reformulated  and  solved  at  each  control  sample  time. 
While  much  of  Equation  (28)  thus  seems  to  be  “wasted”,  it  is  necessary  to  point  out 
that  the  “horizon  of  prediction”  L  has  a  very  strong  effect  on  closed  loop  performance. 

(2)  The  arrays  defining  U*  can  rapidly  become  highly  unwieldy.  This  is  especially 
true  in  most  LSS  applications.  However,  two  mildly  restrictive  assumptions  can  be 
exploited  to  significantly  reduce  required  core: 

(i)  keep  NSM  fixed  in  value  for  each  of  the  L  future  control  settings; 


Q  =  diag(Qu...,QL) 


where  Qi  is  a  pNSM  x  pNSM  diagonal  submatrix  of  Q. 

Assumption  (i)  was  discussed  earlier,  and  probably  has  no  limiting  effect  of  any 
significance.  Assumption  (ii)  is  potentially  restrictive,  but  is  very  common  in  practi¬ 
cal  design  algorithms.  If  these  assumptions  are  allowed,  the  much  lower  dimension, 
symmetric  mL  x  mL  matrix  ( HTQH )  may  be  created  directly  from 

Uhtqh)>A  =  E  e*h*-,+. 

k=maz(ij) 

where  h,  are  defined  in  Equation  (30).  This  can  be  done  offline,  or  whenever  the  system 
model  is  updated.  The  following  computation  must  be  done  every  update  cycle,  as 
the  greatest  core  savings  require  combining  operations  directly  with  the  latest  outputs 
(although  such  an  inefficiency  in  core  could  be  tolerated  if  processing  speed  is  a  greater 
concern):  the  mL  x  1  vector  HTQ(Yd  —  Y-*,)  may  be  created  from  the  L  elements  of  size 
m  x  1 

£  bTkQk+L-t(Yd  -  Yzi)k+L-i  (33) 

t=i 

where  (Yd  —  Yzi)k  axe  the  pL  X  1  partitions  of  the  p  •  NSM  ■  L-dimension  vector  ( Yd  —  Y*) . 
Equation  (33)  is  a  corrected  and  dimensionally  upgraded  (SISO  to  MIMO)  version  of 
Equation  (20)  in  Reid,  et  al.,  (1981).  Thus  only  the  L  h,  in  Equation  (30)  are  needed 
to  generate  U* ,  not  all  of  H. 

Alternate  solution  methods  to  investigate  are  solving  Equation  (31)  via  the  Singular 
Value  Decomposition  (if  memory  is  not  a  problem),  or  solving  Equation  (28)  directly 
by  such  a  process,  as  long  as  certain  relationships  and/or  limits  are  defined  for  control 
values. 

2.3.2  Closed  Loop  Stability  Analysis 

The  optimal  control  sequence  u *{k)  can  be  conveniently  written  as 

u*(k)  =  L0y(k )  -  L,x{k)  +  Les 


where 


(34) 


I'wsnrv  i 


&  V 


tv  V 


Np  is  the  number  of  output  prediction  points,  and  the  L,  are  the  NSM-L  row  submatrices 
of  Lg, 

Lg  =  [L\,  Li, LmM  l) 

The  Li  are  of  size  m  x  m,  and  La  is  itself  defined  from  Equation  (32): 

«*(*)  =  [/m,0,...,0]£7*  4  IJJ' 

=  Iu(HTQH  +  R)-1HTQ{Yd-Y2i) 

=  Lg(Yd  -  Yzt)  (35) 

where  Im  is  the  m  x  n  identity  matrix.  Using  output  relations  defined  earlier,  the  result 
(35)  leads  directly  to  Equation  (34),  for  a  constant  setpoint,  s(k)  =  s,  and  for 

Vd{k)  =  (Jp  -  A *)a  +  A*y(0), 

a  specific  version  of  Equation  (21). 

Due  to  the  form  of  Equation  (34)  we  can  interpret  L0  as  the  “output  feedback”  gain 
matrix,  L,  as  the  “state  feedback”  gain  matrix,  and  Lc  as  the  “set  point  feedback”  gain 
matrix.  Equation  (34)  also  demonstrates  that  the  MPC-optimal  control  law  is  implicitly 
equivalent  to  a  combination  of  a  state-feedback  and  an  output-feedback  control  law,  and 
the  optimal  control  loop  is  equivalent  to  the  set-up  in  Figure  3.  We  note,  however,  that 
the  sequence  u*(k)  is  best  implemented  via  the  approach  developed  in  the  preceding 
subsection,  which  avoids  the  explicit  need  for  full  state  feedback. 

It  is  now  straightforward  to  compute  the  equivalent  closed  loop  matrix  from  the 
block  diagram  in  Figure  3.  Clearly,  the  closed  loop  states  will  evolve  as 

xu(k  +  l)  =  Fxu{k)+Gu'{k) 

=  Fxet(k)  +  G\L0yci[k)  -  L,xet(k)  +  Les ]  (36) 

=  FdXciik)  +  Gets 

where 

Nr 

Fct  =  F+  J2GLi(KC-CFi), 

t=i 

Np 

Gd  =  ZGLiil-ti), 

»= i 

and  s  is  the  command  set  point. 

Again  we  can  interpret  Fet  as  the  equivalent  closed  loop  matrix  and  Gct  as  the 
equivalent  input  matrix.  The  optimal  closed  loop  input  sequence  u*(k)  and  output 
sequence  y(k)  axe  related  to  xct{k)  as 

u*(k)  =  (L0C  -  L,)xel{k)  +  Les  (37) 

y{k)  =  Cxet(k)  (38) 


Summarizing,  the  MPC-optimal  control  sequence  u*(k)  can  be  generated  in  two 


ways: 


mourn 


Figure  3:  An  MPC- Equivalent  Loop. 


(a)  at  each  k,  compute  the  vector  Z(k )  =  (Yi  —  V„)  and  then  u'(k )  =  LaZ(k)-, 
or 

(b)  compute  the  matrices  L0,  L,  and  Lc  or  obtain  the  optimal  control  sequence 
from  the  setup  in  Figure  3  or  simulate  the  system  derived  in  Equations  (36) 
to  (38). 

Method  (a)  is  preferred  for  on-board  implementation:  it  can  run  at  roughly  the 
same  time  as  method  (b),  but  is  more  core-efficient.  Method  (a)  can  be  made  fully  core 
efficient,  resulting  in  very  substantial  core  savings,  but  at  some  cost  in  processing  speed. 
This  tradeoff  awaits  more  detailed  analysis. 

It  has  been  shown  that  the  MPC  control  technique  is  equivalent  to  a  combination 
of  a  constant  gain  state  and  output  feedback  control  law  and  therefore  can  be  analyzed 
in  a  classical  control  framework.  The  stability,  robustness  and  sensitivity  of  the  MPC 
technique  can  be  ascertained  a  priori  in  this  technique  and  the  designer  can  iterate  on 
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various  design  parameters  to  improve  the  performance  of  the  loop. 

The  model  predictive  control  approach  provides  a  transparent  relationship  between 
the  system  performance  and  design  parameters.  The  inherent  robustness  of  MPC  al¬ 
lows  it  to  change  automatically  some  of  its  control  parameters  on-line  when  there  is 
performance  degradation.  MPC  will  efficiently  generate  stabilizing  control  commands 
for  plants  which  are  open-loop  unstable  and/or  nonminimum  phase.  The  latter  is  often 
a  characteristic  of  LSS  structures  with  failed  or  degraded  components.  The  ability  of 
MPC  to  adapt  optimally  in  response  to  system  or  environment  changes  brings  to  the 
LSS  a  very  intelligent  form  of  active  control.  It  is  also  clear  that  passive  or  “stand-alone” 
control  is  achieved  L  oth  by  the  feedback  nature  and  robustness  of  MPC. 

2.4  Simulation 

Simulation  analysis  was  performed  using  linear  discretized  “truth”  and  control  design 
models.  Section  3  details  the  results  obtained  in  Phase  I  of  this  project.  Early  results 
were  obtained  using  a  laboratory-scale  flexible  model  (Tahk  and  Speyer,  1987),  dis¬ 
cretized  to  a  basic  sample  period  of  0.05  seconds.  This  model  is  dynamically  equivalent 
to  the  disk-torsion  bar  problem  considered  by  Cannon  and  Rosenthal  (1984). 

PC-MATLAB  was  the  environment  for  developing  a  “test”  version  of  MPC.  While 
greatly  restricted  by  size  and  run-time  limitations,  this  has  been  an  excellent  envi¬ 
ronment  for  validating  extensions  of  the  algorithm,  general  debugging,  and  analyzing 
dynamic  behavior.  The  CVA  identification  code  was  not  converted  to  PC-MATLAB 
because  (i)  mature  VAX/ VMS  versions  exist,  and  there  is  little  structural  upgrading  of 
the  CVA  algorithm  needed  for  it  to  be  applied  to  LSS  identification;  and,  (ii)  the  algo¬ 
rithm  is  sufficiently  detailed  that  too  much  Phase  I  effort  would  be  needed  to  convert 
to  PC-MATLAB. 

A  VMS  Fortran  version  of  MPC  also  is  being  used  on  larger  models.  It  utilizes 
features  which  provide  speed  at  the  expense  of  core  -  for  example,  NSM  (see  Section 
2.3)  is  replaced  by  an  L-  dimension  vector  nsm,  thereby  allowing  for  a  variable  number 
of  smoothing  points  per  control  update  over  the  “horizon  of  prediction.” 


3  ANALYSIS  OF  RESULTS 


The  scope  and  resources  of  Phase  I  of  this  project  preclude  an  in-depth  investigation 
of  the  full  set  of  control  design  issues.  We  have  nevertheless  been  able  to  initiate  at 
least  a  preliminary  examination  of  most  of  these.  The  results  presented  here  indicate 
clearly  that  the  identification  and  control  methods  offer  feasibility  for  useful  application 
to  flexible  structures,  but  also  that  more  detailed  analysis  on  more  realistic  models  is 
needed  in  order  to  define  properly  the  role  for  LSSICS.  We  propose  that  this  analysis 
be  expanded  in  Phase  II  of  this  project. 

3.1  Algorithm  Modifications 

Section  2  describes  new  results  in  the  area  of  refining  and  upgrading  the  Model  Predic¬ 
tive  Control  algorithm.  The  major  features  of  this  effort  are: 

(i)  a  fully  MIMO  version  of  MPC  has  been  developed  and  successfully  tested 

(ii)  key  errors  in  the  SISO  algorithm  of  Reid,  et  al.,  (1981)  have  been  eliminated 
in  the  upgraded  version 

(iii)  options  have  been  developed  for  the  following: 

-  core  vs.  run  time,  relating  to  exploiting  certain  mildly  restrictive  as¬ 
sumptions  about  the  design  parameters  NSM  and  Q  to  greatly  reduce 
the  core  required  to  generate  the  control  gain  matrix  (cf,  Equation  (32), 

Section  2.3) 

-  formulating  elements  of  H  either  in  terms  of  impulse  response  or  step 
response  functions 

-  mechanization  based  on  input /output  relations,  or  upon  closed  loop 
stability  matrices  (see  (iv)) 

(iv)  for  certain  formulations  of  output  quantities,  analytic  expressions  for  closed 
loop  MPC  have  been  developed  which  allow  for  standard  MIMO  stability  and 
performance  analysis  (See  Section  2.3).  It  is  thus  now  possible  to  observe 
directly  the  effect  of  changes  in  MPC  design  parameters,  and  also  in  the 
control  design  model.  The  stability  analysis  segment  of  MPC  is  meant  to 
be  performed  offline  (though  it  could  become  an  online  part  of  accepting  or 
changing  design  parameters  following  certain  types  of  model  updates).  Thus, 
it  is  not  as  core  efficient  as  the  input/output  version,  but  does  produce  the 
same  response  results. 

It  is  also  very  likely  that  closer  examination  of  the  particular  formulations 
(see,  for  example,  Equations  (34)  and  (35))  will  indicate  a  way  to  develop 
core-efficient  relations  of  the  stability  analysis  code. 

The  last  point  in  item  (iv)  above  also  applies  to  the  core/run-time  tradeoff  concerning 
whether  NSM  is  constant  or  not.  Core  savings  and  algorithmic  simplifications  result 
when  NSM  is  constant;  however,  there  is  definitely  a  way  to  save  most  of  the  core  now 
being  “wasted”  for  the  variable  NSM  version  of  MPC.  To  exploit  this  savings  would 
greatly  add  to  the  matrix  index  bookkeeping  (arrays  would  have  to  be  added  to  the 
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code  just  for  this  purpose)  and,  as  stated  earlier,  the  benefits  to  Phase  I  of  such  an 
undertaking  are  at  best  unclear. 

Some  of  the  core/run-time  issues  are  quantified  in  the  following  analysis  summarized 
in  Table  1.  In  the  example,  we  have  used  NSM  =  4,  L  —  10,  and  m  =  p  =  5.  Distinct 
storage  is  needed  for  all  quantities  identified  in  the  table;  also,  most  “core-efficient” 
quantities  are  submatrices  of  the  “run-time-efficient”  matrixes.  Using  the  above  num¬ 
bers,  representing  a  rather  modest  LSS  design  case,  we  observe  a  clear  difference.  The 
core-efficient  code  requires  13%  of  the  storage  needed  by  the  other  formulation,  al¬ 
though  we  note  again  that  there  is  a  way  to  bring  the  latter  much  closer,  to  perhaps 
about  twice  as  much  required  storage.  At  this  level,  run-time  comparisons  would  be 
quite  meaningful. 


Table  1:  MPC  Software  Implementation 


Issue 

Storage 

Run-Time 

Degree 

VS. 

pr  =  [uf, 

Roughly  equal 

uT  =  K.Auf,...] 

Minor 

Constant 

vs. 

Varying  NSM 

Constant  NSM 

Varying  NSM 

Large,  either  way 

H  using  step  response 

vs. 

Impulse  /  response 

Step  response 

Roughly  equal 

Minor 

Td  Value 

No  effect 

Inversely  Related 

Moderate 

L  Value 

Directly  Related 

Directly  Related 

Moderate 

At  this  point,  run-time  vs.  storage  tradeoffs  can  be  made.  Table  1  compares  storage 
data  for  MPC,  and  similar  data  exists  for  CVA.  More  or  less  core  will  be  required 
depending  on  the  following  issues  summarized  in  Table  2.  Specific  test  scenarios  must 
be  run  to  quantify  the  classifications  cited  in  the  table.  As  stated  earlier,  higher-level 
language  CVA  code  is  by  itself  rather  completely  optimized.  Note  that  it  will  be  possible 
to  do  much  to  offset  the  excessive  storage  advantage  held  by  the  fixed-NSM  option,  if 
sophisticated  matrix  indexing  logic  is  added.  To  see  that  this  is  so,  observe  that  H 
has  L  submatrix  rows,  to  contain  the  “constant  control  submatrices,”  h.  With  varying 
NSM  for  each  of  the  L  future  control  settings,  the  L  submatrix  rows  of  H  each  have 
a  different  “height.”  Thus,  there  are  L(L  —  1)  distinct  h  matrices  (see  also  Equation 


Table  2:  Storage  Comparison  of  Different  MPC  Implementations 


Key  Elements  of  MPC  Gain 
Matrix,  Control  Computation 

Run-Time  Efficient 
(no  Assumptions  on  NSM,  Q) 

Core-Efficient: 

NSM  =  constant,  Q  diagonal 

General 

Example 

General 

Example 

La  =  (fr'QH  +  R)-1  HTQ 

(L  ■  m)  x  (NSM  Lp) 

10,000 

( L  ■  m)(L  ■  m  +  1) 

2,550 

Yd,  Y„ 

2[(NSM  •  L  ■  p)  x  1| 

400 

2[(NSM  •  L  ■  p)  x  Ij 

iTM 

R 

(NSM  ■  Lp)  x  (L  ■  m) 

10,000 

L  x  (NSM  p  m) 

Q 

(NSM  ■  L  ■  p)2 

40,000 

L  x  (NSM  •  p)  x  NSM  ■  p) 

BPlillll 

R 

(L  •  m)  x  (L  x  m) 

2,500 

Lx  m2 

250 

TOTALS 

62,900 

8,200 

(29)): 


[h[‘> 

0 

...  o' 

H  = 

hi11 

h«,” 

: 

: 

0 

.hi" 

...  h[L\ 

(39) 


Each  has  NSM(fc)  mxp  response  matrices  (cf,  Equation  (30)),  for  1  <  k  <  L.  Much 
core  can  be  saved  by  noting  that,  for  given  fc  and  nsm  =  mint(NSM(fc)),  the  first  nsm 
submatrices  in  each  hj^  are  equal,  for  all  j:  that  is  hj^  =  hj^  =  •  •  *,  through  their  first 
nsm  submatrices.  The  issue  to  be  quantified  here  is  whether  the  payoff  in  performance 
merits  such  an  effort. 

All  issues  of  this  nature  will  be  analyzed  and  documented  further  in  Phase  II  in 
detail  sufficient  to  allow  firm  decisions  based  on  operational  priorities. 


3.2  CVA  Identification  Results 

To  demonstrate  the  above  methods,  a  system  dynamically  equivalent  to  the  disk-torsion 
bar  problem  considered  by  Cannon  and  Rosenthal  (1984)  was  used  initially  for  simula¬ 
tion.  This  system  has  a  number  of  similarities  with  large  space  structures  but  is  simple  in 
form  for  simulation  and  analysis.  The  system  consists  of  4  unit  weight  masses,  or  disks, 
arranged  along  a  line  and  connected  by  torsion  bars  modeled  by  spring  and  damper 
elements.  The  system  input  is  a  torsion  on  mass  2  and  the  output  is  the  measured 
rotation  displacement  of  mass  2.  All  of  the  disks  are  free  to  move  about  their  common 
axis  (Figure  4).  The  system  is  an  8  state  system  with  modes  at  0.0  for  the  rigid  body 
mode  and  at  1.4721,  1.9219,  and  2.1598. 

For  demonstrating  CVA,  a  white  noise  of  standard  deviation  of  10  was  used  as  the 
system  input,  and  the  output  was  measured  in  the  presence  of  a  white  measurement 
noise  with  various  standard  deviations.  The  sample  rate  was  0.01  seconds  over  10 
seconds  for  a  total  of  1000  measurements.  The  CVA  was  done  using  40  lags. 
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The  results  of  the  system  identification  are  illustrated  in  Figures  5  and  6.  The 
transfer  function  frequency  response  is  shown  for  the  true  system  and  for  the  identified 
system  for  the  cases  of  white  measurement  noises  of  standard  deviations  0.0,  0.0005, 
0.005,  and  0.05  respectively.  The  system  identification  accuracy  depends  upon  the  mag¬ 
nitude  of  the  measurement  noise.  Note  that  the  identification  is  excellent  for  the  lower 
measurement  noise.  As  the  noise  increases,  the  ability  to  resolve  the  lower  frequencies 
diminishes.  A  plot  with  linear  frequency  would  show  that  the  lower  frequency  peaks  are 
progressively  narrower  and  thus  more  difficult  to  identify  in  the  presence  of  white  noise. 
The  frequency  response  phase  shows  accuracy  comparable  to  the  magnitude  with  high 
accuracy  for  low  measurement  noise. 


Table  3:  AIC  Behavior  with  State  Order 


State 

Order 


Number  of 
Parameters 


Measurement 

0.005 

2504.84 

-4990.22 

-5125.18 

-5119.79 

-5271.53 

-5399.03 

-6226.97 

-6221.00 

-6270.39 

-6265.29 

-6259.64 

-6255.55 

-6249.71 

-6243.97 

-6243.69 

-6237.97 

-6234.42 

-6228.51 

-6222.66 

-6217.65 


Noise  Standard 
0.0005 

2504.83 

-7498.35 

-7802.97 

-8511.80 

-8514.71 

-8830.58 

-8827.75 

-9074.95 

-9601.93 

-9596.16 

-9590.67 

-9586.03 

-9584.38 

-9580.52 

-9574.63 

-9571.49 

-9565.60 

-9559.61 

-9554.36 

-9548.41 


Deviation 

0.0001 

2504.83 

-7219.62 

-7631.33 

-9016.18 

-9034.42 

-9259.73 

-9263.09 

-9487.54 

-10471.00 

-10466.10 

-10461.30 

-10455.60 

-10449.70 

-10443.90 

-10438.00 

-10432.00 

-10426.30 

-10420.40 

-10414.40 

-10408.40 


Three  similar  cases  were  simulated  with  slightly  different  values  for  the  measurement 
noise  variance.  The  identified  order  of  the  system  in  all  of  these  cases  was  order  8  which 
was  in  fact  the  true  system  order.  The  values  of  the  AIC  for  these  three  cases  is  shown 
in  Table  3.  Note  that  the  AIC  falls  sharply  until  order  8  and  then  slowly  increases  at  a 
linear  rate  proportional  to  the  number  of  additional  parameters  needed  for  each  increase 
of  1  in  the  state  order. 


3.2.1  Beam  Model 


A  finite  element  analysis  was  conducted  on  a  homogeneous  beam,  with  both  ends 
unattached,  the  so-called  “Free-Free”  Beam  Model.  This  model  is  potentially  very  much 
more  complicated  than  the  disk-torsion  bar  model  above,  due  both  to  its  dynamics  and 
to  the  ability  to  create  a  very  large  system  order.  We  will  briefly  describe  development 
of  a  5-  segment  version  of  this  beam,  and  then  a  2-segment  version.  The  latter  is  almost 
too  coarse  to  be  useful,  but  it  was  necessary  to  start  at  this  level  due  to  the  strain  on 
computing  resources  which  the  5-segment  version  creates. 


Figure  4.  Laboratory  Four-Disk  System 
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Figure  6.  Transfer  Function  Frequency  Response  Phase 
(  -  true  model;  identified  models:  +  0.05,  •  •  •  0.005, _ 0.0005,  —  0.0  ) 

The  general  equation  for  transverse  deflection  of  a  beam  is  given  by 

-  J^j  [£/(*)  ^1  +  /(*,<)  =  (40) 

where  f{x,t)  is  the  applied  force,  t  is  time,  x  is  axial  direction,  and  y  is  beam  deflection 
from  the  axis.  Applying  the  usual  solution  methods  for  the  free-free  boundary  condition 
case  leads  to  the  following  expressions  for  (symmetric.)  elemental  mass  and  stiffness 
matrices: 


Mi  = 


rrijhj 


K,  =  k: 


156  22 h  54  -Uh 
?  Ah?  13/i  -3  h2 

156  -22 h 
Ah? 

12  6/i  -12  6/i 

Ah2  -6 h  2h? 

12  -6  h 

Ah2 


where 


*  =  EL 

1  h) 

E  is  Young’s  modulus,  I  is  beam  inertia,  and  hj  is  the  segment  length. 


Th  j  finite  element  model  segments  each  have  2  nodes,  and  at  each  node  the  trans¬ 
verse  displacements  are  expressed  as  a  translation  and  a  rotation  (Figure  7a).  The  basic 
5-segment  model  developed  for  LSSICS  analysis  consists  of  equal  length  segments  and 
6  nodes.  The  displacement  vector  is  this  of  dimension  12  (Figure  7b). 

The  dynamic  relation  (39)  then  becomes 


Mu  +  Du  +  Ku  =  / 


(41) 


If  we  define  q  as  the  state  variable  vector  of  modal  deflections,  Equation  (40)  converts 
to 

-M'lD  - M~lK  1  f  M-1 
I  0 


9  = 


q  + 


0 


/ 


(42) 


where 


a 

Q  = 


u 

u 


(43) 


and  D,  M  and  K  are  built  from  the  elemental  matrices  My  and  Kj  in  the  usual  way. 
Except  where  indicated,  we  assume  no  damping  (D  =  0).  Note  that  q  has  24  elements, 
and  that  the  submatrices  in  Equation  (41)  are  12-by-12. 

This  model  is  currently  so  large  for  efficient  analysis  of  the  CVA  identification  algo¬ 
rithm  that  a  smaller,  2-segment  model  has  been  developed.  The  vector  q  in  this  case 
has  dimension  12.  A  2-input,  2-output  system  has  been  defined  by  sensing  translation 
and  rotation  deflections  only  at  node  3  (u5  and  u6),  and  placing  actuators  at  node  2  (uj 
and  u4). 
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Figure  7.  Displacement  Variable  Definition  for  5-Segment  Free-  Free  Beam 


3.2.2  Results 

The  first  set  of  results  for  the  application  of  CVA  identification  to  flexible  beam  dynamic 
systems  is  presented  for  the  case  of  the  symmetric  2-segment  beam,  derived  from  the  left 
two  segments  of  the  beam  pictured  hi  Figure  7b.  The  dynamic  model  to  be  identified 
is  given  by  Equation  (41),  and  the  ime  2-input,  2-output  formulation  as  defined  above 
is  also  used  here. 

The  results  are  shown  in  Figure  8  in  the  form  of  transfer  function  magnitude,  mag¬ 
nitude  squared  (“power”),  and  phase  plots  vs.  frequency.  In  this  formulation,  there 


are  four  total  transfer  functions  between  all  combinations  ot  the  2  inputs  and  2  outputs. 
They  are  numbered  like  the  four  elements  of  a  2-by-2  matrix,  with  element  (i,j)  con¬ 
necting  input  j  to  output  i.  The  plots  in  Figure  8  are  representative,  and  do  not,  for 
example,  show  all  of  che  phase  plots. 

The  true  transfer  function  has  two  poles  influenced  by  the  translational  force  input 
u3,  and  another  two  poles  influenced  by  the  torque  input  u4.  The  rigid  body  modes  of 
course  occur  at  zero  frequency. 

Analysis  of  the  data  shows  that  CVA  does  a  good  job  of  locating  the  frequencies 
of  vibration,  but  does  less  well  in  isolating  totally  the  response  of  the  force  input, 
(.,  1)  to  the  three  nodes  of  the  true  system.  That  is,  the  poles  associated  with  the 
torque  input  (.,2),  appear  both  in  the  torque  response,  (2,.),  and  the  force  response, 
(1,.).  This  effect  is  best  seen  in  the  transfer  function  magnitude  plots  of  Figure  8;  the 
transfer  function  squared  plots  tend  to  obscure  this  effect,  and  hence  highlight  the  good 
agreement  between  the  true  and  estimated  responses.  The  representative  phase  plots 
in  Figure  8  are  more  accurate  than  appears,  due  to  the  “wrap-around”  effects  at  ±n. 

Table  4  shows  that  the  identified  basic  system  dynamic  modes  match  those  of  the 
truth  system  rather  well  (the  discrete  system  eigenvalues  are  shown  here;  continuous 
system  comparisons  would  be  equivalent,  given  knowledge  of  the  sampling  interval). 

The  coupling  between  the  torque  inputs  and  force  outputs  mentioned  above  is  felt 
to  be  due  to  the  unrealistically  high  degree  of  symmetry  of  the  original  two  element 
beam  model,  with  the  control  inputs  located  exactly  in  the  center  of  the  beam.  We  thus 
analyzed  the  affects  of  asymmetry  by  developing  a  model  whose  central  node  is  at  the 
60%  point  along  the  beam.  Figure  9  shows  that  the  asymmetry  greatly  enhances  the 
accuracy  of  the  identification.  Plots  (a)  and  (b)  of  this  figure  even  show  fine  agreement  at 
the  smaller  frequencies  in  a  log  plot,  which  is  very  difficult  to  achieve,  since  the  statistical 
“strength”  of  an  identification  is  spread  more  or  less  evenly  over  the  frequency,  and  not 
its  log. 

Plots  (c)  through  (r)  in  Figure  9  plot  the  (a)  and  (b)  results  on  the  same  frequency 
scale  used  in  Figure  8.  Here,  there  is  no  need  to  use  the  square  of  the  transfer  function 
to  highlight  the  similarities.  All  of  the  transfer  functions  agree  very  well  with  the  true 
system  transfer  functions,  both  in  magnitude  and  phase,  as  the  other  plots  in  Figure 
9  show.  Table  5  shows  that  the  modes  for  the  asymmetric  two  segment  beam  are  also 
very  accurately  identified.  The  symmetric  beam  identification  does  at  least  as  well  in 
this  aspect,  as  discussed  above,  which  is  the  significant  information  needed  for  adaptive 
control. 

Identification  was  also  done  on  a  “51%”  beam,  with  results  being  clearly  better  than 
the  symmetric  case,  but  not  quite  as  good  as  the  “60%”  case.  The  conclusion  is  that 
some  asymmetry  is  beneficial  to  some  aspects  of  the  overall  identification  process,  but 
is  not  critical  to  achieving  good  real  time  control. 

Finally,  effects  of  measurement  noise  on  the  identification  are  analyzed.  The  “60%” 
beam  was  used  to  generate  the  results  shown  in  Figure  10.  Plot  (a)  represents  the 
case  cov(R)  =  diag(5.SE  —  06,1.2 E  —  05),  where  R  is  the  symmetric  measurement 
noise  covariance  matrix  (off-diagonal  terms  not  given).  Plot  (b)  represents  the  case 
cov(R)  =  di ag (0.0020, 0.0035).  Each  of  these  plots  compares  directly  with  Figure  9a, 
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Table  4:  System  Eigenvalues,  2  Segment  Beam 


tt 

v 

v  r.' 
1 

K 

V  i'v 
O  V 


J:  i 


fe  $ 


ing,  has  a  magnitude  and  phase  spectrum  as  shown  in  Figure  11.  This  corresponds  to 
a  SISO  system,  with  collocated  sensor  and  actuator  at  mass  2;  the  other  masses  are 
totally  free  to  move. 

The  key  MPC  design  parameters  are:  Td,  system  quantization  interval  (seconds),  r, 
the  vector  of  reference  trajectory  time  constants  (seconds),  error  and  control  weights  Q 
and  R,  smoothing  points  per  control  cycle,  NSM,  and  number  of  control  updates  in  the 
horizon  of  prediction,  L. 

Figure  12  shows  the  result  of  the  MPC-controiied  disk-torsion  bar  system  when 
mass  2  is  commanded  to  slew  2  “units”  from  a  quiescent  state.  In  this  as  in  all  time 
response  plots,  the  abscissa  shows  time  in  seconds.  For  relatively  little  design  effort,  a 
very  effective  result  has  been  obtained.  The  mass  is  controlled  to  its  desired  position  to 
within  1%  error,  after  reaching  this  level  in  under  a  second.  Note  that  the  uncontrolled 
masses  still  exhibit  very  lightly  damped,  oscillatory  responses,  which  require  control 
response  due  to  coupling.  This  “tight”  solution  results  from  zero  control  weighting 
(R  =  0).  A  tradeoff  is  set  up  between  tight  response  and  control  energy  required  for 
same,  as  seen  in  Figure  13,  which  shows  results  for  R  —  0.1  and  1.0.  Note  the  expected 
result  of  large  demands  on  an  unweighted  control  input.  This  study  highlights  the 
ability  of  MPC  to  control  specific  structural  nodes. 

Following  the  closed  loop  stability  analysis  approach  outlined  in  Section  2.3,  the 
appropriate  Fci  and  Gd  matrices  were  created.  Figure  14  shows  how  the  closed  loop 
transfer  function  for  R  =  0  is  superior  to  both  the  R  =  1.0  and  open  loop  case  (latter 
in  Figure  11). 

An  offline  design  procedure  may  clearly  be  developed  based  on  Fel  and  Gci  using 
standard  analysis  tools.  Since  these  matrices  are  explicit  functions  of  the  MPC  design 
parameters,  the  latter  may  be  varied  until  curves  such  as  Figure  14a  result.  This 
analysis,  however,  is  very  burdensome  on  the  core  available  to  a  PC  or  even  a  /i-VAX; 
we  propose  in  Phase  II  to  use  the  stability  matrices  as  both  design  and  analysis  tools 
on  larger  systems,  when  Grumman’s  facilities  are  available. 

This  investigation  will  also  focus  on  ways  to  generate  Fci  and  Gci  more  efficiently. 
For  example,  if  Q  is  purely  diagonal  -  the  only  case  so  far  analyzed  -  there  is  great 
potential  for  very  large  savings  in  core,  since  ( p  •  NSM  ■  L)(p  ■  NSM  •  L  —  1)  locations 
become  superfluous. 

A  final  note:  it  is  apparent  that  the  closed  loop  stability  matrices  Fci  and  Gei  can  be 
used  to  simulate  the  MPC  closed  loop  response.  Doing  so  produces  a  response  which, 
as  expected,  matches  perfectly  the  Figure  12  response. 

3.3.1  Free-Free  Beam  Design 

Analysis  began  first  with  the  2-  segment  finite  element  model  of  the  free-free  beam 
described  in  Section  3.2.  The  same  2-by-2  input/output  structure  defined  in  that  section 
is  used  here  as  well.  This  is  an  exceptionally  difficult  system  to  control,  for  these 
reasons: 

(i)  sensors  and  controls  are  not  collocated 


(ii)  usually,  every  displacement  state  is  controlled  directly;  in  this  case,  only  2 
controls  are  applied  to  6  displacement  states 

(iii)  the  basic  model  has  no  damping  (damping  is  supplied  in  some  cases,  however) 

Item  (iii)  is  emphasized  because  it  makes  the  open  loop  system  technically  unstable. 
The  “difficulties”  imposed  are  partly  due  to  a  desire  to  severely  test  MPC,  but  also  to 
core  limitations. 

Figure  15  shows  the  open  loop  (uncontrolled)  response  of  the  undamped  2-segment 
model  to  a  deflected  initial  state,  no  controls.  This  same  initial  state  is  used  for  all 
examples  in  this  section.  This  is  a  regulator,  as  opposed  to  tracking  problem.  The  four 
“non-output”  responses  represent  the  unmeasured  nodal  deflection  variables,  consisting 
of  two  translation  and  two  rotation  variables.  None  of  the  6  deflection  rate  variables 
are  considered  for  measurement  or  control  in  this  analysis,  although  this  is  clearly  a  key 
item  for  future  work. 

One  of  the  best  MPC  designs  to  date,  which  is  a  clear  improvement  over  the  uncon¬ 
trolled  case,  is  shown  in  Figure  16.  The  fact  that  control  activity  has  lessened  consid¬ 
erably  after  5  seconds  while  not  completely  eliminating  response  oscillations  -  though 
they  appear  to  be  slowly  damping  out  -  indicates  that  this  particular  sensor /actuator 
structure  is  not  optimally  suited  for  control.  Notwithstanding  this  observation,  and  in 
view  of  the  three  design  difficulties  imposed  on  MPC  by  this  particular  problem,  MPC 
clearly  improves  performance.  For  this  case,  R  =  0,  NSM=12,  and  L  =  8,  and  the 
weight  Q  is  a  pure  diagonal  whose  elements  monotonically  increase  from  upper  left  to 
lower  right  on  the  diagonal.  When  the  control  prediction  window  L  is  increased  to  12 
in  this  problem,  the  results  are  better  still,  as  shown  in  Figure  17. 

More  than  the  other  MPC  design  parameters,  it  seems  that  NSM,  L,  and  R  have 
the  most  direct  effect  on  a  good  MPC  design.  Halving  Td  to  0.025  actually  had  a 
destabilizing  effect  on  performance,  as  seen  in  Figure  18.  The  parameter  t  was  increased 
by  a  factor  of  10  to  1.50  seconds  in  another  study,  but  this  had  very  little  effect  on 
response.  Also,  as  pointed  out  by  Reid  et  al.,  (1981),  setting  R  =  0  usually  results  in 
tighter  response.  This  is  of  course  preferred  as  long  as  there  is  sufficient  control  energy. 
We  noted  the  same  effect  in  the  disk-torsion  bar  problem  discussed  above. 

System  damping  remains  an  effective  means  of  reducing  the  system’s  internal  energy 
and  with  it,  the  stabilizing  control  effort  required.  Figure  19  shows  a  sequence  of  MPC 
solutions  for  a  system  in  which  damping  increases  from  none  (a)  to  less  than  about 
0.1%  (b),  to  about  4%  (c).  Solution  (a)  represents  the  first  3  seconds  from  the  more 
complete  solution  in  Figure  16.  (Economy  of  simulation  time  is  very  important:  each 
solution  in  Figure  19  takes  over  an  hour  to  generate  on  a  PC;  the  longer  one  in  Figure 
16  takes  about  6  hours.)  Note  that  with  even  the  addition  of  the  most  mild  damping, 
the  control  effort  decreases  much  more  than  the  output  response  differs.  It  is  not  until 
rather  significant  damping  is  added  (Figure  19c)  do  the  outputs  change  markedly.  The 
open  loop  4%  damping  case  is  shown  in  Figure  21  for  comparison  with  Figure  19c. 
See  also  the  plots  for  the  unmeasured  (“non-output”)  variables  in  Figure  20.  Damping 
greatly  affects  their  amplitude.  This  analysis  supports  Reid’s  rule  of  thumb  that  NSM 
aind  L  eaich  be  at  least  of  order  n. 
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The  next  series  of  simulations  deal  also  with  the  two-segment  beam  model,  but  with¬ 
out  damping.  The  above  analysis  was  deliberately  set  up  to  be  very  difficult  on  a  control 
algorithm;  there  were  much  fewer  inputs  and  outputs  configured  than  there  are  system 
modes,  and  the  basic  translational  and  rotational  displacement  and  velocity  coordinates 
which  have  been  used  generate  a  high  degree  of  numerically  inefficient  coupling.  It  was 
felt  that  transforming  the  physical  coordinates  to  the  "eigen”,  or  modal,  coordinates 
would  provide  MPC  a  better  environment  for  a  good  design. 

The  transformation  is  performed  on  the  dynamic  system  (41)  (with  D  =  0  here), 
and  converted  into  the  expression 


f 


r?  +  At?  =  YTf  (44) 

where  Y  =  Q~lU ,  Q  is  the  Cholesky  factor  of  the  mass  matrix,  M  =  QTQ,  and  U  is 
the  unitary  matrix  resulting  from  applying  the  Schur  decomposition  to  K  —  Q~T KQ~l , 
i.e.,  UAUT  =  K.  Initially,  the  output  equation  was  retained,  as  follows: 


y  =  Cq  =  CYfj  (45) 

where  f?  is  comprised  of  r?  and  r?  in  the  same  manner  as  q  in  Equation  (43). 

The  Figure  17  case  was  rerun  using  the  modal  formulation  just  described.  The 
results,  shown  in  Figure  22,  clearly  show  significant  improvement  resulting  from  the 
modal  formulation.  We  emphasize  again  that  this  is  a  system  with  6  dynamic  modes, 
subject  only  to  2  remote  actuators  and  2  colocated  sensors. 

Somewhat  surprisingly,  when  the  actuators  are  also  colocated  at  one  end  of  the  beam, 
a  poorer  design  results  (Figure  23).  This  is  likely  due  to  a  decrease  in  controllability 
resulting  from  the  distance  of  these  actuators  from  the  other  nodal  points.  Finally, 
when  four  actuators  are  devoted  to  the  control  of  the  end-beam  displacement  variable, 
making  a  4  input,  1  output  configuration,  Figure  24,  the  result  is  even  better  them  the 
case  shown  in  Figure  22  (the  four  inputs  here  are  at  the  3,  4,  5,  and  6  positions  in 
Figure  7,  and  the  single  output  is  beam  displacement  at  station  5).  This  result  matches 
expectations. 

An  MPC  design  has  been  initiated  for  the  24  order  five-segment  free-free  beam  model 
described  in  Section  3.2.  At  the  time  of  this  report,  this  analysis  is  not  yet  complete. 
System  order  has  initially  forced  smaller  values  for  NSM  (4)  and  L( 7)  than  desirable.  To 
offset  this  somewhat,  the  4  sensors  and  actuators  in  this  configuration  were  collocated  as 
shown  in  Figure  25.  There  is  no  damping  here.  The  objective  again  is  to  regulate  each 
output  to  zero  from  a  nonzero  initial  energy  state.  The  two  responses  shown  in  Figure 
26  are  for  the  case  Qi  ~  10(*2)  and  R  =  I.  This  value  of  R  is  driven  less  by  proper 
application  of  design  principles  learnt  thus  far  than  by  accommodating  the  algorithmic 
numerical  deficiencies  of  the  n-V AX  resident  Fortran  code.  At  this  time,  the  more 
numerically  stable  and  core  efficient  MPC  code  is  on  the  PC,  written  in  PC-MATLAB. 
The  Fortran  code  is  currently  being  upgraded. 

This  examples  does  point  out  however,  the  potential  for  MPC  to  control  moderately 
large,  lightly  damped  systems. 
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3.3.2  Robustness  Analysis 

Robustness  is  the  ability  of  a  control  system  to  perform  as  optimally  as  possible  when 
unmodeled  changes  occur  in  the  system,  or  when  the  dynamic  environment  changes  in 
an  unpredictable  way.  MPC  is  a  control  method  known  in  past  applications  to  be  rather 
robust.  A  quantified  analysis  of  the  robustness  properties  of  MPC  will  be  conducted 
during  Phase  II. 

A  short  series  of  simulations  designed  to  investigate  the  robustness  of  MPC  to  spe¬ 
cific  flexible  structures.  In  one  case,  using  the  disk-torsion  bar  model  of  Cannon  and 
Rosenthal  (1984),  a  “failure”  was  generated  at  4  seconds  into  a  simulation  which  in 
every  other  respect  matched  the  “nominal”  case  discussed  above  (see  Figure  12).  For 
this  case,  the  viscous  damping  between  mass  1  and  mass  2  was  halved.  Recall  that 
the  sensor  and  actuator  are  collocated  at  mass  2,  so  that  its  motion  is  to  be  controlled 
directly.  Figure  27  shows  that  MPC  successfully  offsets  the  sudden  removal  of  damping 
within  2  cycles.  Note  also  that  both  nominal  and  “failed”  period  mass  2  trajectories  are 
controlled  to  the  commanded  level  to  a  very  high  tolerance  when  the  control  weight  R 
is  set  to  zero  (Figure  27d).  If  R  is  nonzero,  Figure  28,  that  is,  not  well-matched  to  the 
problem,  the  sudden  system  change  can  result  in  an  unacceptably  oscillatory  response. 

Robustness  of  the  free-free  beam  model  was  analyzed  by  increasing  by  50%  the 
stiffness  of  the  left  segment  (in  the  2-segment  case)  at  5.0  seconds.  Figure  29  shows 
that  MPC  again  is  able  to  reestablish  control. 

In  both  cases,  the  system  or  “truth”  model  was  changed  but  not  the  control  design 
model.  MPC  is  thus  ignorant  of  the  new  situation,  but  is  able  to  perform  well.  It  is 
useful  to  note  that  onboard  identification  can  reestablish  for  MPC  the  proper  control 
design  model.  This  is  a  very  important  feature  of  LSSICS;  robustness  of  the  “new” 
configuration  is  clearly  enhanced  by  accurate  knowledge  of  the  system  dynamics. 

As  this  phase  of  the  project  was  ending,  some  work  has  begun  on  the  effect  of  actua¬ 
tor  loss  (output  force/torque  =  0  following  damage)  on  system  performance.  The  MPC 
results  discussed  above  were  deliberately  generated  using  fewer  control  inputs  than  sys¬ 
tem  modes.  A  practical  design  would  be  more  conservative  in  this  area,  but  this  analysis 
does  demonstrate  the  ability  of  MPC  to  stabilize  and  provide  robust  performance  for 
“under-controlled”  systems.  Some  comparison  was  done  with  Linear  Quadratic  designs, 
using  full  state  feedback.  Although  the  results  of  the  comparison  are  preliminary,  the 
LQ  designs  seem  to  be  more  sensitive  in  performance  to  the  loss  of  an  actuator.  For  a 
baseline  LQ  design  utilizing  the  full  six  actuators,  failing  only  one  of  these  at  a  time, 
in  the  two  segment  free  beam  case,  actually  causes  instability  in  four  of  the  six  cases. 
The  MPC  results  discussed  above  are  clearly  more  robust,  although  it  is  recognized  that 
stable  LQ  designs  are  possible  using  only  two  controls.  Further  effort  is  needed  in  this 
analysis  before  firm  conclusions  can  be  made. 

3.4  More  New  Results 

This  section  presents  results  generated  too  late  for  proper  location  in  this  report. 

A  control  design  scheme  called  Independent  Modal  MPC  has  been  developed  to 
enhance  the  effects  and  benefits  of  modal  control  which  were  described  above.  IMMPC 


operates  on  the  notion  that  if  a  system  is  controllable,  it  can  be  controlled  by  controlling 
each  mode  separately.  An  MPC  design  is  performed  on  each  mode  sequentially,  and 
then  the  resulting  controls  are  combined  into  a  final  command  input.  Figure  30  shows 
the  result  when  this  process  was  applied  once,  to  the  highest  frequency  mode  of  the  two 
segment  beam.  The  figure  shows  the  difference  in  the  response  of  this  mode  when  it 
is  uncontrolled  (solid  line)  and  controlled  (hatched  line).  Note  that  better  closed  loop 
damping  seems  achievable  using  this  approach. 

A  modal  design  using  IMMPC  was  next  performed  configuring  sensors  of  displace¬ 
ment  (or  rotation)  and  velocity  at  each  coordinate  of  the  2  segment  beam  model.  Mea¬ 
surements  and  control  inputs  are  synchronous  at  100  Hz,  and  the  MPC  weighting  ma¬ 
trices,  Q  and  R,  were  adjusted  appropriately  for  this  case.  The  results  (Figures  31 
and  32)  show  that  MPC  seems  to  work  very  well  when  there  is  sufficient  monitoring 
of  all  dynamic  modes.  Figure  31  presents  the  uncontrolled  (hatched  line  here,  opposite 
from  Figure  30!)  and  controlled  (solid)  modal  space  responses,  and  Figure  32  the  cor¬ 
responding  physical  space  responses.  The  length  of  time  in  these  simulations  is  very 
brief,  0.2sec.,  but  appears  adequate  to  show  the  desired  decay  of  the  responses  to  zero. 

It  is  of  course  necessary  to  carry  this  analysis  further  by  designing  reduced  order 
observers,  and  also  by  investigating  effects  of  sensor  failures.  However,  it  does  seem 
that  the  basic  questions  about  the  applicability  of  MPC  to  LSS  control  are  favorably 
answered. 
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(m)  Estimated  Transfer  Function  -  (1,2)  Element 


(n)  True  Transfer  Function  -  (1,2)  Element 
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(o)  Estimated  Transfer  Function  -  (1,2)  Element 


(p)  True  Transfer  Function  -  (1,2)  Element 
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Estimated  Transfer  Function  -  (2,2)  Element 


(r)  True  IVansfer  Function  -  (2,2)  Element 


Figure  8.  CVA  Analysis  of  12th  Order  Symmetric  Free  Beam 
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(g)  Estimated  Transfer  Function  —  (1,2)  Element 


(h)  True  Transfer  Function  —  (1*2)  Element 


(m)  Estimated  Transfer  Function  -  (2,1 )  Element 


(n)  True  Transfer  Function  -  (2,1)  Element 


d)  Estimated  Transfer  Function  -  (2,2)  Element 
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(q)  Estimated  Transfer  Function  -  (2,2)  Element 


(p)  True  Transfer  Function  -  (2,2)  Element 
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(r)  True  Transfer  Function  -  (2,2)  Element 


Figure  9.  CVA  Analysis  of  12th  Order  Asymmetric  Free  Beam 
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(a)  Measurement  Noise  Covariance  =  diag{5.8E  -  06, 1.2E  -  05) 


ESTIMATED  TRANSFER  FUNCTION  -  (1.  1)  ELEMENT 
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(b)  Measurement  Noise  Covariance  =  dtag(0.0020, 0.0035) 


Figure  10.  Effect  of  Measurement  Noise  on  CVA  Identification 
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(a)  Beam  Responses,  NSM,  L=12 
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(b)  Beam  Inputs,  NSM,  L=12 
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(c)  Beam  Non-Output  Responses,  NSM,  L=12 


Figure  17.  Two  Segment  Free  Beam  Closed  Loop  Response  Using  MPC;  L  =  12 
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(a)  Beam  Response,  Modal  4  Input,  4  Output  Formulation 


(c)  Force  and  Torque  Inputs  at  Output  Station 


(b)  Force  and  Torque  Inputs  at  Beam  Center 
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(d)  Non-Output  Responses 


(e)  First  Z  Non-Output  Responses  (Detail) 


Figure  24.  Two  Segment  Free  Beam  response,  Using  4  Inputs  and  1  Output 
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27.  Robustness  Analysis  of  Disk-Torsion  Bar  System;  Damping  Reduced  by 

50%  at  4  Seconds 


mode  I  response 


mode  5  response 


displacement  displacement  displacement 


JS 

s* 

& 


SM 


t 


0  0.05  O.i  0.15  0.2 


0  0.05  O.i  0.15  0.2 


Figure  32.  Controlled  vs.  Uncontrolled  Physical  Responses  Using  Independent  Modal 

MPC  Design  Method 
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4  SUMMARY  AND  CONCLUSIONS 

4.1  Summary  of  Results 

The  results  presented  above  indicate  clearly  that  CVA  identification  and  Model  Pre¬ 
dictive  Control  are  capable  of  providing  robust  control  of  flexible  structures.  This 
conclusion  is  the  result  of  applying  CVA  and  MPC  to  three  flexible  structure  models: 
an  8th  order  disk-torsion  bar  system  with  one  input  and  one  output,  collocated;  a  12th 
order  free-free  beam  with  2  inputs  and  2  outputs;  and  a  24th  order  free-free  beam  with  4 
inputs  and  4  outputs.  The  sensors  and  actuators  were  not  collocated  in  the  12th  order, 
2-segment  beam  problem.  Both  the  smaller  and  larger  beam  models  were  analyzed  with 
basically  zero  damping,  and  stabilizing  control  was  achieved.  For  a  sufficient  level  of 
input  noise  and  also  for  the  proper  level  of  lags,  and  some  process  and/or  measurement 
noise,  CVA  does  well  in  system  identification,  although  more  analysis  and  tuning  of  the 
larger  beam  case  remains  to  be  done. 

The  analysis  to  date  has  established  that,  for  the  identification  and  control  of  flexible 
structures,  the  horizon  of  prediction  must  be  adequately  large.  Specifically,  L  should 
be  at  least  the  order  of  the  system,  and  so  should  the  number  of  output  smoothing 
points  when  the  system  has  minimal  or  no  damping.  Both  parameters  may  be  reduced 
somewhat  for  a  structure  with  as  little  as  2%  damping.  Also,  it  was  determined  that 
the  system  sample  period  is  a  sensitive  parameter,  as  is  the  direct  control  weight  R. 
Assuming  adequate  control  resources,  best  performance  is  achieved  for  R  =  0.  The 
dimension  size  of  the  output  error  weight  Q  is  very  important,  but  performance  thus 
far  is  rather  insensitive  to  changes  in  the  values  of  Q  elements. 

Sensor/ actuator  location  is  important.  Deliberately  difficult  configurations  were 
analyzed  in  the  beam  problem,  and  analysis  points  to  the  need  for  careful  consideration 
of  the  configuration.  MPC  stability  analysis  can  assist  in  this  task. 

Inherent  robustness  of  LSSICS  using  MPC  was  demonstrated,  and  the  potential  for 
its  enhancement  using  CVA  identification  discussed. 

Further  robustifying  at  the  algorithm  level  may  result  from  performing  the  design 
in  terms  of  modal  coordinates,  which  decouple  the  input/output  structure,  or  in  terms 
of  the  statistically  significant  dynamic  states  as  determined  offline  using  CVA. 

The  analysis  done  in  Phase  I  has  confirmed  many  suppositions  but  has  naturally  un¬ 
covered  or  left  unresolved  other  important  issues,  many  of  which  we  propose  to  address 
in  Phase  II. 

4.2  Conclusions  and  Recommendations 

The  work  performed  under  this  phase  of  the  project  indicates  clearly  that  the  combina¬ 
tion  of  CVA  identification  and  adaptive  control  using  MPC  offer  potential  for  being  a 
key  part  of  robust  LSS  control  synthesis  methods.  There  are  further  details  to  resolve, 
but  results  obtained  to  date  indicate  that  it  is  worthwhile  to  pursue  their  resolution. 
With  regard  to  the  effectiveness  of  CVA,  the  conclusion  is  that  some  asymmetry  is  bene¬ 
ficial  to  some  aspects  of  the  overall  identification  process,  but  is  not  critical  to  achieving 
good  red  time  control. 
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Recommendations  fall  into  the  area  of  performing  detailed  analysis  and  comparison 
of  this  approach  with  other  methods.  Metrics  related  to  LSS  mission  requirements 
and  performance  specifications  should  be  developed  and  incorporated  into  the  design 
process.  More  work  is  needed  in  fully  integrating  the  CVA  and  MPC  algorithms.  Also, 
further  analysis  of  the  sensor  and  actuator  configuration  and  its  effect  both  on  a  useful 
control  design  and  on  system  performance  should  be  done.  Location  of  actuators  and 
sensors  can  be  combined  into  the  MPC  design  process  using  nonlinear  programming 
methods. 

The  benefits  and  costs  of  decentralized  control  of  flexible  orbiting  structures  should 
also  be  analyzed,  particularly  from  the  standpoint  of  CVA  and  MPC  applications.  Fi¬ 
nally,  means  should  be  found  to  incorporate  well-understood  metrics  such  as  control¬ 
lability  and  observability  into  the  selection  of  key  MPC  design  parameters,  e.g.,  the 
weighting  matrix  elements. 
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LSS  Analysis  Software 

Grumman  Corporation  has  several  programs  for  large  scale  control-system  synthe¬ 
sis  and  analysis  as  well  as  large-scale  structural  analysis,  which  assist  in  comparison 
analyses. 

Optimal  Control  Codes.  Grumman’s  extensive  library  of  efficient  and  accurate 
optimal  control  and  estimation  computer  codes  enable  the  user  to  design,  analyze, 
simulate,  and  evaluate  large-scale  problems  in  a  relatively  short  time.  The  codes  reside 
on  permanent  libraries  on  both  the  batch  and  time  sharing  systems  so  that  they  are 
readily  accessible.  Two  major  control  analysis  systems  are  used. 

DIGISYN.  DIGISYN  is  a  design-oriented  program  which  performs  direct  analog 
or  digital  designs  and  analyzes  the  results  via  covariance  analysis,  linear  time  history 
simulations,  frequency  response  analysis,  etc.  It  allows  the  user  to  perform  partial 
state  feedback,  full  state  feedback,  and  optimal  estimator  designs.  DIGISYN  is  used 
extensively  in  aircraft  and  spacecraft  design  efforts.  Grumman  developed  this  program 
using  corporate  funds  and  considers  it  proprietary. 

CASCADE.  A  second  set  of  codes,  this  is  a  library  of  proprietary  Grumman  pro¬ 
grams,  to  be  used  extensively  in  this  study,  will  be  invaluable  for  large  space  structures 
research.  Some  of  the  unique  features  of  this  package  are: 

•  A  built  in  capability  to  do  normal  mode  analysis  for  systems  with  many  degrees 
of  freedom  (from  several  hundred  to  a  thousand) 

•  A  capability  to  specify  the  sample  time  required  to  control  the  system  based  on 
the  overall  uncertainty  which  will  be  encountered  compared  with  tolerable  control 
system  deviations 

•  A  capability  to  design  a  continuous  or  digital  control  systems  which  minimizes  the 
same  criteria 

•  A  numerical  approach  which  does  not  require  integration  to  compute  equivalent 
discrete  systems  and  performance  indices  (all  of  these  depend  on  the  use  of  a 
generalized  eigenvalue  eigenvector  program) 

•  The  computation  of  the  optimal  full  state  control  gain  or  Kalman  filter  gain,  in 
steady  state,  for  continuous  or  discrete  designs  which  uses  the  “Potter”  approach 
(again,  this  is  an  eigenvalue-  eigenvector  approach) 

•  The  computation  of  the  “reduced  state”  feedback  gain  in  continuous  or  discrete 
time  which  solves  for  the  optimal  gain  using  a  Davidon  minimization  of  the  cost 
equation 

•  A  capability  to  design  Luenberger  observers  of  any  specified  dimension  and  the 
specification  of  the  minimal  observer  dimension 

•  A  general  controllability-observability  analysis  which  permits  sensor/actuator  place¬ 
ment  so  that  control  authority  and  measurement  information  are  optimally  dis¬ 
tributed  among  the  various  states  being  controlled 

•  A  general  maximum  likelihood  or  recursive  filtering  approach  to  parameter  identi¬ 
fication.  This  program  is  independent  of  the  control  system  design  code  that  does 


r„ 

*  * 


,N 

*  * 
.V 


rj 


all  of  the  previous  analysis/synthesis. 

CASCADE  enables  the  user  to  develop  the  design  model  directly  on  the  interactive 
terminal.  He  specifies  the  nearest  mode  in  the  finite  element  model  that  the  sensor 
is  mounted,  the  type  of  sensor  (linear  or  rotation,  absolute  or  differential),  where  it  is 
located,  and  whether  it  is  a  force  or  torque  device.  It  applies  forces,  the  direction  that 
the  force  is  applied.  For  external  disturbances,  the  user  specifies  the  locations  of  the 
disturbance  (distributed  disturbances  appear  at  many  or  possibly  ail  nodes  on  the  finite 
element  model).  The  matrices  A,B,C  and  D  in  the  model: 

x  =  Ax  +  Bu  (46) 

y  =  Cx  +  Du  (47) 

are  then  determined  by  CASCADE  by  reading  the  appropriate  data  from  a  NASTRAN 
structural-analysis  output  file. 

Analysis  and  Simulation,  both  linear  and  nonlinear  simulation  packages  of  vary¬ 
ing  degrees  of  complexity  have  been  developed  at  Grumman.  The  linear  package  allows 
time  history  generation  by  two  methods. 

1.  Propagation  of  difference  equation 
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where  $  is  the  state  transition  matrix,  and  the  matrices  A  and  B  are  defined  above. 

2.  Integration  of  the  differential  equation  between  sample  times. 
PROTOBLOCK.  PROTOBLOCK  is  a  Grumman-developed  program  that  enables 
the  user  to  create  a  control-system  block  diagram  at  a  computer-graphics  workstation. 
Each  block  can  represent  an  elementary  control  element  or  a  higher-order  transfer  func¬ 
tion.  Once  the  block  diagram  has  been  created,  PROTOBLOCK  can  automatically 
generate  the  state  equations  as  well  as  time-  history  graphs  of  the  response.  The  block 
diagram  can  be  easily  modified  and  the  simulation  can  be  quickly  repeated  to  assess  the 
impact  of  design  changes. 

SATSIM  (Satellite  Simulation).  This  Grumman-developed  program  enables  an 
extremely  detailed  simulation  of  a  satellite  attitude  control  system.  To  facilitate  the 
computation  procedure,  small  relative  structural  motions  are  treated  in  a  linear  way 
while  nonlinear  simulation  is  employed  for  large  rigid-body  motions.  Features  of  this 
program  include: 

•  Built  in  capability  to  handle  detailed  spacecraft  configurations 

•  Ability  to  simulate  many  control  system  configurations 

•  Detailed  models  of  actuators  and  sensors 

•  Nonlinear  modeling 

•  Accurate  computation  of  expected  disturbances 

•  An  automatic  NASTRAN  interface 

SATSIM  is  in  extensive  current  use  to  support  work  for  both  military  and  civilian 
spacecraft  design  efforts. 


Structural- Analysis  Programs.  Grumman  has  many  structural-  analysis  com¬ 
puter  programs;  among  them  are  DISCOS,  SPACE14,  NASTRAN,  ASTRAL-COMAP 
and  PLANS.  DISCOS  (Dynamics  Interaction  Simulation  of  Controls  and  Structure)  is 
designed  to  analyze  the  stability  of  controlled  spacecraft  which  consist  of  coupled  flex¬ 
ible  bodies,  SPACE14  is  a  Grumman-generated  program  for  the  analysis  of  controlled 
rotating  satellites.  A  class  of  deployment  problems,  on-board  mass  motion  (e.g.,  crew 
members)  and  fluid  motion  can  be  simulated,  and  environment  disturbances  (gravity 
gradient,  solar  radiation  pressure,  and  aerodynamic  loads)  are  automatically  generated. 
Both  the  MacNeal-Schwendler  version  of  NASTRAN  and  Grumman ’s  structural  ana¬ 
lyzer,  COMAP  ASTRAL  are  used  extensively  for  creating  and  analyzing  large-  order 
finite-element  models.  PLANS  (Plastic  Analysis  of  Structures)  is  a  Grumman  devel¬ 
oped  program  used  for  the  analysis  of  structures  which  can  have  plastic  deformations. 
It  has  been  used  extensively  in  computerized  crash  analysis  of  automotive  and  aircraft 
structures. 
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